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Optimal Resource Allocation for Two Processes* 


By C. COURCOUBETIS*t and P. VARAIYA* 
(Manuscript received May 9, 1984) 


Two processes compete for access to n resources. A scheduling policy 
allocates the resource when the processes request it simultaneously. The 
objective is to minimize the average value of a state-dependent cost. The 
optimal policy can be calculated explicitly for the case of one resource. In the 
general case n > 1, an adaptive scheduling algorithm is proposed. The algo- 
rithm measures average transition times of the system and converges to the 
optimal policy. 


I. INTRODUCTION 


Two processes share n resources. A process operates in one of these 
states: thinking, requesting, or holding a given resource. The thinking 
and resource holding times are geometric, with means depending on 
the process and the resource. 

Time for the system is discrete. A decision has to be made when the 
two processes simultaneously request the same resource k, in which 
case a scheduling policy assigns the resource to process 1 with proba- 
bility u, € [0, 1]. A resource cannot stay idle when there is a process 
requesting it. The problem is to choose u* € [0, 1]", which minimizes 
the average value of a cost depending on the state of the system. 

A simple example of such a system is the case of two processes 
sharing the same broadcast facility. Such a process goes through the 
following phases. It thinks for an arbitrary amount of time (any 
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activity other than broadcasting), requests the resource (broadcast 
channel), uses the resource for an arbitrary time (broadcasts), and 
then resumes thinking. If the resource was busy upon request, the 
process has to wait until the resource is released. If two processes 
request the resource simultaneously (during the same slot), the sched- 
uler of the resource will decide, according to the scheduling policy of 
the system, which process will proceed first. 

In Section III the problem is solved explicitly for the case n = 1 and 
when the cost is taken to be resource idle time. It is shown that for 
large ratios of average thinking times, u* does not depend on the 
holding times and contrasts with the familiar “cu” results, where the 
average waiting cost is minimized when priority is given to processes 
with large waiting cost and small holding times (e.g., Refs. 1 through 
6). It turns out that in that case the optimal policy consists in serving 
the process with the least thinking time first. 

This is a new result for nonpreemptive systems. Similar results 
appeared in the literature concerning n processes sharing a single 
resource in a preemptive-resume way. For such systems, the “least 
thinking time first” scheduling rule has been proved optimal (see Refs. 
7 and 8). For example, in Ref. 8 the author considers a “mirror image” 
problem, where the goal is to minimize the utilization (repair time) of 
_ the scarce resource (repairman). In this case the policy of serving the 
process with longest thinking time first is optimal. Because preemptive 
policies were considered and a single resource was shared, results are 
not applicable to the problem considered here. 

In Section IV we consider the case of n resources and an arbitrary 
state-dependent cost. Since Section III suggests that an analytic 
expression for the optimal policy u* as a function of the parameters 
of the processes would be extremely hard or even impossible to get, 
we proceed with an adaptive algorithm to compute u*. This algorithm 
converges in a finite number of steps to u* by using measurements of 
the average transition times of the system. The analysis is based on 
the results contained in Refs. 9 through 15, which are reviewed in 
Section 4.1. In Section 4.4 the structure of the problems for which the 
results apply is generalized. We permit a process’to consist of an 
arbitrary set of thinking states, an arbitrary set of resource holding 
states for each resource in the system (one set of states per resource), 
and a number of resource request states (one state per resource). 
Finally, in Section V we present some open problems. 


“Il MODEL DESCRIPTION 


In this section we give a formal description of the model for n = 1. 
The state of process 1 (i = 1, 2) at time t (t = 0) is X}. If X; = 0, 
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then process i is thinking. If Xi! = 1 [resp. 2], then process i is 
requesting [resp. holding] the resource. See Fig. 1. 

The resource can accommodate only one process at any time. When 
the two processes are simultaneously requesting the resource, then the 
resource is assigned with probability u € [0, 1] to process 1. 

Denote by p; the probability that process 1(1 = 1, 2) will complete 
its thinking time in a given time unit. Similarly define q; for completion 
of service (holding). From this definition, one finds that the process 
X, = (Xi, X?), t = 0, 1, ---, is a Markov chain with transition prob- 
ability matrix P = {P(x, y)|x, y © X = {0, 1, 2}*} defined as follows: 


P(00, 10) = pi(1 — pe), P(00, 01) = (1 — pi)po, P(00, 11) = pips; 
P(OL, 12) = pi, P(10, 21) = pe; 

P(01, 02) = (1 — p,), P(10, 20) = (1 — pg); 

P(11, 21) = u, P(11, 12) =1-u; 

P(02, 12) = pi(1 — gz), P(02, 00) = (1 — pidge; 

P(20, 21) = (1 — qi)p2, P(20, 00) = qi(1 — pr); 

P(12, 20) = qo, P(21, 02) = qu. (1) 


The diagonal elements of P are defined so that the rows sum to one. 
For any choice of u, X; is ergodic. 

Define as z,, the invariant measure of X;, and let k: X — R be the 
state-dependent cost vector. Then the expected cost-per-unit time is 








SCHEDULER-CONTROLLED 
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7 


Gj 
REQUEST STATE 


THINKING STATE 


Fig. 1—Process i (i = 1, 2). 
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T 
Iu) = lim TE Y RX) = Y malx dela). (2) 


t=0 xEX 


Our goal is to find u* that minimizes (2). 


Il. MAXIMIZING RESOURCE UTILIZATION 


Throughout this section we will assume k(00) = 1, k(y) = 0, y 4 00. 
Let z be the invariant probability measure associated with P. Then 
the expected idle time of the resource is 


J(u)= Y mali). (3) 


i,7E{0,1} 
Direct computation of z,, shows that (3) takes the following form: 


Au+B 
Cu + D’ 


where A, B, C, D are functions of {p, pe, gi, G2}. From (4) it follows 
that 


J(u) = (4) 


sign (224) = sign(AD — BC). (5) 
In what follows, we will choose u* € [0, 1], which minimizes (4). From 
(5) it follows that u* is an extreme point of [0, 1]. We state now the 
key results, which are proved in the Appendix. 


Fact: 
AD — BC = PipeqigsF iF, 


where F; = F;(pi, P2; 41, q2); l = 1, 2. 

Theorem 1 implies that the sign of AD — BC is determined by the 
sign of F;. 
Theorem 1: 


F,20 forall py, pe, dQ, ge € [0, 1]. 


We state now the main theorem. 
Theorem 2: 1. The curve {(p1, P2) € [0, 1]?| Fi(pi, pe, G1, G2) = 0} lies 
between the lines p; = pz and p2 = p;/2 for qi S qez, and between p, = Pe 
and p, = P2/2 for q2 < qu. 
2. p; = 2pe implies F, = 0, and pz = 2p; implies F, < 0, for all qiq2 © 
[0, 1]. 
The following corollary follows now from Theorem 2 and (5): 
Corollary 2.1: Let u* be the value of the control minimizing (4). Then 
the following hold: 

1. py = 2p2 implies u* = 1, and pez = 2p, implies u* = 0; 
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0 1 Dy 


Fig. 2—Optimal policy u* as a function of the parameter py, po. 


2. If pi = Po, then gz S q, implies u* = 0, else u* 

3. If q1 = qa, then p; = P2 implies u* = 1, else u* 

Theorem 2 and Corollary 2.1 are illustrated in Fig. 2, where F, = 0 
is considered for some fixed values of qi, q2, 2 < qi. The surprising 
result here is that there exist “safe” regions of the parameters, where 
u* does not depend on the service times of the processes. 


= 0. 


IV. ESTIMATION OF u* 


In this section we give a method for calculating u* for an arbitrary 
cost vector k. The algorithm we propose requires no a priori knowledge 
of the parameters of the processes and is performed adaptively to the 
system. It starts by applying some arbitrary policy u°, and then by 
monitoring the evolution of the system and estimating certain tran- 
sition times, it updates the policy used until u* is reached; this is 
always achieved in a finite number of steps, depending on the number 
of the shared resources. Since the algorithm is adaptive, it could be 
effectively used in systems with slowly varying parameters. This, 
together with the simplicity of the computation involved, makes this 
approach an interesting alternative to a direct calculation of u* by 
using standard policy iteration algorithms (see Ref. 15). In Section 4.1 
some general results are stated. We discuss in Section 4.2 the case of 
two processes and one resource, and in Section 4.3 the case of two 
processes and n resources. 


4.1 Some results in Markov decision theory 


Propositions 1, 2, and 3 consist of results already known and are 
stated without proof in order to make the present work self-contained. 
The key propositions used to calculate u* are Propositions 4, 5, and 
6, which consist of new results; their proofs are included in the 
Appendix. 
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Let X; be a Markov chain on {1, ---, s}. If X, = 1, then any control | 
u € U(it) may be used, where U(i) is a compact set. A stationary 
strategy is any element u = (u(1), ---, u(s)) € U = U(1) X --- xX 
U(s). Let P(u) be the s X s transition probability matrix describing 
X;, and assume that X, consists of a single ergodic class and that P(u) 
is continuous on U. Let Q(u) := P(u) — I, and let k := (kR(1), ---, 
k(s))’ be the cost vector, not depending on u. The cost to be minimized 
is J(u) = lim, 1/(T' + 1)E bo k(X;). The following two propositions 
give the optimality conditions for u. 
Proposition 1: (see, e.g., Ref. 14, Lemma 3.1) For u € U consider the s 
linear equations in the s + 1 variables y € R,c € R*, 


yl = Q(u)c + k, (6) 


where L = (1, ---, 1)’. Then: 

1. If (y, c) is a solution, then y = J(u). 

2. If (y, c) is a solution, then so is (y, c + 61) for all 6. 

3. A solution always exists and is almost unique in the sense of (2). 

Let H(c, u) = Q(u)c + k, H = (AM, ---, H,)’. Note that H,(c, u) = 
H;(c, u(t)) depends only on u(i). Let h(c) = (hi, ---, hs), hi(c) = min 
{Hi(c, v) |v € U(i)}. 
Proposition 2: (see, e.g., Theorems 3.1 and 3.2) The control u is optimal 
(minimizing J(u)) iff there exist (y, c) such that yl = h(c) = 
H(c, u). 

The following propositions will be used in the part dealing with the 
estimation of u*. 

Consider the equation 


V, = (I — aP(u)) 1k. (7) 


Then V,(i) is the expected discounted cost with the discount factor a, 
a € [0, 1], starting from state i, ie., Va(t) = ES 20 atk(X;)}. Let 
Z? := (I — aP(u))™. 

Proposition 3: (see, e.g., Ref. 9, Chapter 3) Z° is the fundamental 
matrix of the absorbing chain X? defined on {0, 1, ---, s} with the 
(s + 1) X (s + 1) transition probability matrix P, such that 


P,(i, J) = aP(i,j), 47 € {1,---, 8}, 
P.i,0)=1-—a for 1#0, P,(0, 0) =1. 
Let N%{[j] be the expected number of visits to state 7 starting from 1 of 


X; in a geometrically distributed interval of time with mean (1 — a)’. 
Then z% = N$[J]. 
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Proposition 4: 
lim[V.(i) — Va(R)] = c(i) — c(k), 
a—1 


where c satisfies (6). 
Proposition 5: Let T;; be the expected time before the first visit to state 
j by X; starting from 1. Then 
Ty = Ty See 
“HY ig RAFI 
lim[N 71] — Ni{jl] = 
if kR=j Al. 


oJ 


i 
The next proposition is key, since it relates the dual variable c to 


transition times of the chain X;, which can be estimated quite readily. 
The proposition can be readily obtained from Propositions 3, 4, 5, and 


eq. (7). 
Proposition 6: 


==> a ty k(j) + R(i) — k(R). 


j=1 


4.2 Optimal resolution of conflict for two processes and one resource 


In this section we provide the following results. In Theorem 3 we 
characterize the value of u*. In Theorem 4 we relate the value of u* 
with the sign of a quantity that can be estimated using Proposition 6 
from the behavior of the system. We finally propose Algorithm 1, 
which uses these results to adaptively calculate u*. 

Theorem 3: u* can always be restricted to the set {0, 1}. 

Proof: Let (y*, c*) be the variables in (6) corresponding to the optimal 
policy u*. Then by Proposition 2 and the fact that u enters only in 
the row of P corresponding to the state (11), 


min{—c*(11) — (1 — u)c*(12) 
+ uc*(21) + k(11)|u € U} 
min{[c*(12) — c*(11) + k(11)] 
+ u[c*(21) — c*(12)]|u € U} (8) 


and the minimum is achieved at u*. Let A = c*(21) — c*(12). Then 
A > 0 implies u* = 0, A < 0 implies u* = 1, and A = 0 implies that 
any u € (0, 1] willdo. O 

Theorem 4: Let (vy, c) be a solution to (6). Then pienicte) — c(12)] = 
const. for all u € U. 


* 


ne 


RESOURCE ALLOCATION 7 


Proof: Suppose there exists a Up such that c'°(21) — ¢°(12) = 0. Then 
y, c“ are optimal dual variables since the optimality condition (8) is 
trivially satisfied. Then every u € U will satisfy the optimality condi- 
tions with c*(21) — c*(21) = c”(21) — c“(12) = 0; hence u is also 
optimal. From this and Proposition 4.1 it follows that c”(21) — 
c41(12) = c“”(21) — c%(12) = 0, for all u, € U. This, with the next 
Fact, proves Theorem 4. 
Fact: The vector c of (6) can be chosen to be a continuous function of 
u. 

This can be proved as follows. Let 7, be the invariant probability 
measure of P(u). Since y = 7k, by using (6) we get 


TRL - k= Q(ude. (9) 


Since z,, is continuous in u (by the ergodicity of P(u)), and Q(u) is of 
rank s — 1 for all u € [0, 1], then there is always a solution ¢ = (¢(1), 
-++,é(s — 1), 0) of (9) continuous in u, and by the “almost” uniqueness 
of c the Fact follows. O 

Corollary 4.1: Let (y, c) be a solution to (6). Then c(21) — c(12) > 0 
implies u* = 0; c(21) — c(12) < 0 implies u* = 1; and if c(21) — 
c(12) = 0, any u* will do. 

Corollary 4.1 suggests the following algorithm to estimate u* adap- 
tively. 
Algorithm 1: 

1. Start the system with an arbitrary u. 

2. Use Proposition 6 to estimate the sign of c(21) — c(12). 

3. Use Corollary 4.1 to choose u*. 

We will conclude this section with an example. Consider the case of 
minimizing the probability of conflict, i.e., the probability that both 
processes request simultaneously. In this case, k(t) = 0 for i ¥ 11. It 
follows that 


Ti211 ~ T2111 


c(21) — c(12) = R(11) Ta 


and 
sign[c(21) — ¢(12)] = sign[Ti211 — Toiail- 


An intuitive justification for the above equation is the following. 
Minimizing P(11) is equivalent to maximizing T};11, which is equiva- 
lent to giving priority to the process so that the busy period of the 
system corresponding to that process being serviced first is maximized. 
It is easy to see by using renewal arguments that this is equivalent to 
choosing the largest among the 7211, T21,11. 
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4.3 The case of n resources 


In this section we will generalize some of the previous results. Each 

process has a thinking state 0 as before, and pair (1, 2,) of request 
and resource holding states for every resource k. Hence there are n 
conflict states, and the control u is the n-tuple (uw, ---, u,), where wu, 
is the probability of assigning resource k to process 1 in the conflict 
state (1,, 1,). See Fig. 3. Again, we are concerned with the estimation 
of u*. To simplify notation, we will prove the theorems for n = 2; the 
same proofs hold for any n larger than 2. We follow the sequence of 
the previous section. 
Theorem 5: u* can always be selected from {(0, 0), (0, 1), (1, 0), (1, 1)}. 
Proof: Assume that each process consists of the states (0, 11, 2, 
1,, 22), where states 1,, 2; are associated with resource 1, and ly, 2) 
with resource 2. The proof is similar to the proof of Theorem 3, since 
the optimality conditions are 


* 


Y 


min{[e*(112;) — e*(1111) + R(111,)] 

+ uy[e*(231,) — c*(1122)] Jur, € U4}, (10a) 
min{[c*(1222) — c*(1el2) + R(1e12)] 

+ ugfc*(2212) — c*(1222)][ue € U2}. O (10b) 


The following algorithm is a policy iteration algorithm that starts 
with an arbitrary initial choice of u and in a finite number of repeti- 
tions converges to u*. Step 2 corresponds to the value determination 
operation, and Step 3 to the policy improvement operation. The novel 
feature of this algorithm is the simple and adaptive execution of the 
value determination operation by using Proposition 6. The ergodicity 
of the system ensures that every improvement of the policy corre- 
sponds to a strict decrease in cost (see Refs. 15 and 16). Define A, := 





Fig. 3—Process 1 (i = 1, 2), n resources. 
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c(21) — c(12) and Ag := c(2’1’) — c(1’2’). (Observe that A;, As are the 
coefficients of u,, U2 in eq. [10].) 
Algorithm 2: 

1. Start with some arbitrary u* = (uj, w4). 

2. Using Proposition 6, estimate the sign of A? and A, where A? = 
A;(u*), l= 1, 2. 

3. Choose in an extreme way the u/'t!, i = 1, 2, in order to decrease 
each Afu‘*! separately (i.e., if A;(u*) > 0, choose u#*! = 0, ---). If 
u**! % u*, go to Step 2 by using u = u**). If u**! = y*, then u* = u*. 

The convergence follows, since every choice of a different u* corre- 
sponds to a strict decrease of the cost (by ergodicity), and there is a 
finite number of choices for u* (finite extreme points in U). The 
essential difference between the cases n = 1 andn > 1 is that Theorem 
4 cannot be generalized for n > 1. This prevents us from inventing a 
“one step” algorithm for the estimation of u*. 

Note that although not addressed in this paper, the problem of 
estimating the 7;,’s is of great importance for the algorithm. Substan- 
tial errors in the estimation procedure could lead to a wrong choice of 
u* by the algorithm. 


4.4 A generalization 


One can notice that throughout Sections 4.2 and 4.3, the only parts 
of the structure of P used were the rows corresponding to conflict 
states. This leads to a generalization of the form of the processes. A 
process can consist of an arbitrary set Sp of thinking states and a pair 
(1,, S%) for each resource k, where S$ is an arbitrary set of service 
states. The only constraint on the transition diagram of the process is 
that there is a unique state in S§ to which a process k can transit from 
state 1,. 


V. CONCLUSIONS 


One open problem is the relation between the distributions of service 
and thinking times, and the safe regions of Section III. In other words, 
is there a general rule suggesting that for large ratios of thinking time 
of the processes, the choice of u* is independent of service times? This 
would be nice, since no further calculations are needed to obtain the 
optimal u*. 

Another open problem is the generalization to m processes sharing 
n resources. A realistic model for this situation would suggest a 
decentralized information structure for each resource scheduler. By 
this we mean that the scheduler of resource r should base its decision 
on “local” information only, i.e., the state of resource r and the identity 
of processes requesting resource r. If one uses Markov Decision Theory 
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as was done here, the optimal decision of scheduler r concerning 
conflict between processes 1, j would depend on “global” information 
about the system, i.e., the state of all other processes even if they are 
not involved in this particular conflict. This is unsatisfactory. Instead 
one would like to obtain a rule for making local decisions that are 
optimal “in the average”, by “smoothing out” what happens in the rest 
of the system. 

As a final open question, we state the generalization of Section III. 
One should be able to prove the existence of the safe regions of 
Theorem 2 without explicitly calculating the invariant probability 
measure. 
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APPENDIX 
Proof of Theorem 1: After calculating, we have 


AD — BC = pip.qig?F iF», 
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where 
Fy = pipeq2 — Piped — Pid2 — PiP3Q. — Pip3M + pin 

+ pip — p3 — pip, — po+pitpi (11) 
F, = qige[pip3 + pipe — PiP2 — Pr — P2 + 1] 

+ (po — Pip2)g2 + (Pi — Pipe)”. (12) 
Let L(pi, p2) = pipe + pip — Pipe — Pi — P2 + 1. Then by (12) 
PF. = qig2L + (p2 ~ PiP2)g2 + (Pi — PiP2)q1- 

Consider G;(qi) := Fo( pi, D2, G1, G2). Then 


Gi(qi) = qilgoL + (pi — pipe)] + (De — Pipe)qo. 


Since G;(q;) is linear in qi, to prove G, = 0 for all p;, pz, qi, G2 € [0, 1], 
it is enough to show G,(0) = 0 and G,(1) = 0. But G,(0) = (pe — pipe)qe 
= 0; hence we only have to prove G,(1) = 0. Let G2(qz2) := G,(1). Using 
a similar argument, since G2 is linear in gq, and G,(0) = pi — 
D2Pp2 2 0, we only have to show G,(1) = 0. But this holds since 
G2(1) = pipe( pi + Pp2 — 2) — pipe + 1 = 0, as one can readily check. O 
Proof of Theorem 2: Consider the function G,(p2) = Fi(p1, pe, qi, G2) 
and q; < q2. Then proving the theorem is equivalent to proving the 
following: 

1. Gi(p2) has a unique root p? in the interval [p;/2, pi]. 

2. Gi(p2) has no root in the intervals [0, p:/2], [ pi, 1). 

3. D2 < p2 implies Gi(p2) = 0, and ps = p$ implies G,( ps) < 0. 
One can now prove 1 through 3 since 


Gi(1) < 0, Gi(pi) < 0, Gi(pi/2) = 0, Gi(0) = 0, p$p2p3 > 1, 


where p3, p3, p3 are the roots of G,(p2) = 0, and F;(a, 6, c, d) 
—F,(b, a, d, c), as one can easily check. 


Proof of Proposition 4: Define é := (c(1), ---, c(s — 1), 0), and Q := 
(P — I). Then 
yl =Qc+k has a unique solution (vy, ¢). (13) 


Let ya := V,. — V.(s)L. Then by using (7) and subtracting, we get 
(1 — a) V,(s)L = (aP — Dya + R. (14) 


By multiplying (7) by the invariant measure 7, we get z[(1 — a) 
V, — k] = 0, which, together with the ergodicity assumption (ie., all 
components of 7 are >0) and the fact that V, = 0, implies that no 
component of V, tends to © as a > 1. Let y,. := (1 — a) V.(s). By the 
previous argument it follows that there is a converging subsequence 
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of yz as a — 1, and for this subsequence let y* := limg_,.’va. By using 
(14) we find 


Vol = (aP — I)y, + R. (15) 


Consider (15) as a — 1 along the above subsequence. Since it has a 
unique solution y, for every a, it follows that there is a unique y* = 
lim,41¥a, and (13) implies that y = y* and ¢ = y*. The proof now 
follows by observing that V,(i) — V.(k) = y(t) — yo(k), and c(t) — 
c(k) = é(t) — (Rk) by Proposition 1. O 

Proof of Proposition 5: Let X; be the chain under consideration, X; € 
1 = {1, ---, s}, and define Y;, ¢ = 0,1, ---, Y; © Mo := {0, 1} such 
that P[ Yin. = 1] Y; = 0) =1—a, P[Yn1=1| Y: = 1) = 1,1. is an 
absorbing chain with one being the absorbing state. We also define 
the following: 


Nj] := Number of visits by X; to state j starting from 
t before absorption, 


T, := Min{n = 0|X, = J}, 
T, := Min{n = 0| Y, = 1}, 
z := Number of visits to state j before absorption. O 


Lemma: N?{j] = Eyj,olz|P:i[T2 > Ti). 
Proof: Start the system (X;, Y;) from state (1, 0) and count visits to 
state j before absorption. One can always start counting from time T 
on, since no visit to state j occurs before T. The count z will not be 
identically zero only when (X;, Y;) = (Jj, 0), and this occurs with 
probability PT: > T,]. By the strong Markov property, one can 
always restart the system from state (Jj, 0), and the result follows since 
the new expected count will be E;;, 0) [2]. 

Note that N3[j] = Ey, o[z] since P;[T2 > T,] = 1. 

Let 6} := P{X? =] for some t > 0}. Then by using the previous 
lemma, it follows that 


Ni[j] = NfLi\8§, for i Fj, 
and by using similar renewal arguments, 
NiLi] = 1 + BENS; 


hence 


Nj[j] = 





1 
65. 
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Since 6% = E[a™], by using dominated convergence we obtain 


0 - 
a, bo = ElTia™] 
and 
0 
lim — 
sn Oa 


Therefore if 147 4 k, then 
lim [N2Lj] — NeLi]] = lim [(65 - 68) NFU 


BG = Ty. 


= lim By — Bi = By 

a1 1 = By; 
a3 Ti = Tj 
Ty 


by de l’Hospital’s rule. The proof for k =j is similar. O 
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An Infant Mortality and Long-Term Failure Rate 
Model for Electronic Equipment 
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This paper describes the reliability model used by system designers at 
AT&T Bell Laboratories to predict component and equipment reliability. A 
decreasing-failure-rate Weibull model describes the high incidence of early- 
life failures, or infant mortality. This is combined with the constant-failure- 
rate (exponential) model traditionally and widely used for the long term. 
Formal modeling of both early-life and long-term reliability is needed to 
manage the development and manufacture of reliable products. The effects of 
temperature and electrical stress on failure rate are taken into account. A 
model for the effect of integrated circuit dynamic burn-in on reliability is also 
described. 


I. INTRODUCTION 


Reliability describes the ability of a system to continue to perform 
its required function to the satisfaction of the user. Predicting the 
reliability of a new electronic system is an important part of the 
system design process. If the design will not meet reliability objectives, 
it must be improved by using more reliable components, adding system 
redundancy (tolerance to component failures), or performing burn-in 
or other screening. 

The ability to satisfy customers is the most critical factor for a 
viable product. But reliability has other economic impacts as well. 
Repair costs, both during a warranty period and beyond, depend on 
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the reliability. Reliability also affects the numbers of spares needed in 
the field to accommodate expected repair needs. 

System reliability depends on the reliability of its components in 
their application environment. Unfortunately, it is impossible to pre- 
dict the operating life for any individual electronic component. It is, 
however, possible to treat large populations of such components sta- 
tistically with acceptable results. For example, the number of failures 
among a large population of components and each component’s prob- 
ability of survival or failure can be estimated. The statistical behavior 
of the components then determines the statistical behavior of the 
entire system. In this way the reliability of an electronic system can 
be estimated. | 

This paper describes the techniques and models currently used by 
system designers at AT&T Bell Laboratories to predict the failure 
probabilities of electronic components. There are other widely used 
reliability prediction models. Most are based on a simple constant- 
failure-rate model. One such widely used method is described in Ref. 
1, MIL-HDBK-217D. For many purposes, these models are inade- 
quate. Early-life reliability predictions provide information needed to 
balance requirements, design, screening, and field support activities 
for commercial products. By not reflecting the relatively high failure 
rates associated with early equipment life, the constant-failure-rate 
models do not provide the information needed to manage these early- 
life reliability issues. 

The model described here is more realistic. It has elements in 
common with the conceptual “bathtub curve” reliability model. Under 
this model, a relatively large number of defects can be expected early 
in equipment life. This is referred to as “infant mortality”. The | 
likelihood of failure then falls dramatically to a low, constant level 
called long-term reliability. This low failure rate is the behavior 
expected of mature products. Eventually, components can degrade and 
the incidence of failure increases during “wear-out”. The model pre- 
sented here reflects both infant mortality and long-term behavior. 
Integrated circuit reliability dominates the reliability of modern elec- 
tronic equipment. Because properly designed and manufactured silicon 
integrated circuits do not experience “wear-out” behavior, it is not 
modeled. 

To go along with the model, we also have tabulated, elsewhere, 
reliability estimates for a wide variety of components. These estimates 
are based on AT&T Bell Laboratories data whenever possible. Other- 
wise, estimates are obtained from MIL-HDBK-217D, the de facto 
industry standard. 

As we already mentioned, predictions of reliability are useful in 
estimating its impact on both customer satisfaction and economic 
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viability. Both the early-life and long-term aspects of reliability are 
important and need to be addressed. The need to address long-term 
reliability is well known. In addition, formal modeling of early-life 
reliability provides information essential to managing the design, 
manufacture, reliability testing, and screening programs needed to 
assure that initial product reliability will satisfy customers. 


Il. INFANT MORTALITY (SHORT-TERM RELIABILITY) 


Infant mortality is characterized by an initially high, but rapidly 
decreasing, failure rate. The early failures come from a small fraction 
of the components considered to be weak. These weak units contain 
defects (usually manufacturing defects) that are not immediately fatal 
but that will cause failure in a relatively short time. Examples of these 
defects are poor internal electrical connections, the presence of con- 
taminants, and insulating layers that are too thin. 

Failures due to infant mortality can appear in two different ways. 
In one, failures occur during operation after some time. These are 
called “device operating failures” or DOFs. The failures are time 
dependent. The infant mortality part of the reliability model describes 
their occurrence. 

In addition to the DOFs, initial failures are found at various first 
tests, including first circuit pack tests, first system test, or when the 
system is first tested after shipment to the field. These failures are 
called “dead on arrivals” or DOAs. They cannot be related to operating 
time. A component can test as satisfactory, be assembled into equip- 
ment, and then fail to work. No operation has occurred. Instead, these 
failures may be thought of as event dependent rather than time 
dependent. Somehow, handling during equipment manufacture has 
induced failure of the weak component. DOAs are not reflected in the 
reliability model. Although their existence is well recognized, we do 
not know how to quantitatively predict their occurrence. 

During infant mortality, components exhibit a “high incidence of 
failure”, high relative to later life (the long term). This should be kept 
in perspective. Only a very small fraction of components actually fail 
(typically much less than one percent). See the discussion in Section 
VII on calculating numbers of failures for more details. 


Ill. LONG-TERM RELIABILITY 


Infant mortality failures are mostly caused by defects. Even in the 
long term, some failures continue to appear due to manufacturing 
defects; however, other failures occur due to more fundamental com- 
ponent properties. The important aspect, though, is that the failure 
rate is low and relatively constant in the long term. This is the behavior 
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observed in large populations of mature components. Failures occur 
at a fairly constant rate within the entire population; therefore, it can 
be treated as a homogeneous population of components having con- 
stant failure rates. 


IV. RELIABILITY DEFINITIONS 


Before looking at the specific reliability model, we should review 
some basic reliability definitions. Reliability models are based on the 
probabilities of survival or failure of a component or system. These 
probabilities can be described by one of several common functions.” 
Assume that a component starts to operate at time t = 0. Then F(t) 
represents the probability that the component fails at or before time 
t. This is called a cumulative distribution function and it has the 
properties 


F(t) = 0 for t <0, 
0<F(t)< F(t’) for Ostst’, 
F(t) 1 for t— 0, 


The reliability function, R(t), gives the probability of surviving past 
time t. It is related to F(t): 


R(t) =1 — F(t). 


This function is the source of the usual definition of reliability as “the 
probability of surviving”. The derivative of F(t) is a probability density 
function represented by f(t): 


d 
m ft) = 7 FO 


F(t) = { f(x)dx. 


In practice, the instantaneous failure rate or hazard rate \(t) is 
often more useful than the functions just mentioned. From this point 
on, failure rate will mean instantaneous failure rate: 


AY. _ © intro, 


MO RGN a 


which implies that 
RW) =e { Madde 
The failure rate \(¢) of a unit has the following interpretation: If the 
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unit has survived until ¢, the probability of failing in a small time 
interval At at t is At A(t). 


V. COMPONENT FAILURE RATE MODEL 


The reliability model described here applies to individual compo- 
nents. We use the failure rate function to describe the reliability 
model, since it is the most convenient. 

The basic reliability model consists of two parts, shown by the heavy 
lines in Fig. 1. During the infant mortality period, the failure rate is 
described by a two-parameter Weibull model. The Weibull failure rate® 
can be expressed as 


A(t) = Aut. 


This distribution appears as a straight line when plotted on logarithmic 
scales, as in Fig. 1. The slope of the line is —a and the intercept at 
t = 1 hours is ),. The failure rate is initially high but decreases rapidly. 

Beyond 10,000 hours the model assumes that the failure rate is 
constant. The exponential model, which implies a constant failure 
rate, is used. The long-term failure rate is simply 


A(t) = AL: 


Defining the Weibull-to-exponential switchover point to be at ex- 
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Fig. 1—A component failure rate model (shown by solid lines) combining a Weibull 


model and an exponential distribution. The failure rates shown here are typical, but are 
not intended to correspond to any particular component. 
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actly 10,000 hours is arbitrary but reasonable. Weibull behavior has 
been observed to persist for at least a year (8760 hours). At that point, 
the failure rate is changing very slowly. Therefore, any time somewhat 
greater than one year could have been chosen as the switchover point 
with little impact on the modeled failure rates; 10,000 hours was a 
conveniently round number. 

There are two distinct sources of information on component relia- 
bility: direct monitoring of performance in the factory or field, and 
accelerated life tests. Factory or field data give a real measure of 
component reliability in the short term (a year or so), and some data 
are also available for the long term. Accelerated life tests provide 
information about reliability expected in the very long term (tens to 
hundreds of years). 

Reliability studies seldom continue for more than two or three years. 
Therefore, primarily infant mortality is observed. Plotting the loga- 
rithm of the observed failure rate versus the logarithm of operating 
time usually gives a straight line. The straight line means that a 
Weibull distribution describes the failure rate behavior well, as we see 
in Fig. 2. We use the Weibull model to describe infant mortality 
because of such data. 
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_ Fig. 2—Infant mortality removal rate of beam-lead sealed-junction T*L integrated 
circuits (small-scale integration/medium-scale integration) (see Ref. 4). 
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As just mentioned, we can directly measure the time dependence of 
infant mortality failures. However, measuring the lifetime distribution 
of the main population in the long term is impractical using normal 
operating conditions, since so few components fail. Therefore, we use 
accelerated test conditions to estimate the lifetime distribution and, 
hence, the failure rate, at normal use conditions. The test results imply 
that only a very small fraction of semiconductor components will fail 
at normal use conditions within a forty-year service life. This is a very 
important result. It means that wear-out will not occur during the 
service life. (In fact, wearout of semiconductor components should 
never occur.) Figure 3 illustrates this point. For semiconductors, a 
lognormal model describes the main lifetime distribution in accelerated 
tests.” Two such lognormal failure rate curves, extrapolated to normal 
use, are shown in Fig. 3. These examples represent values in the range 
usually observed. 

Accelerated testing does not define the long-term failure rate very 
precisely during the service life. Accelerated life conditions are far 
removed from normal use conditions; therefore, a long extrapolation 
is required to estimate real field performance. This is inherently a 
tricky business. Furthermore, in accelerated life tests, the sample size 
is generally small. In such cases the tests cannot accurately show the 
distribution of the first few percent of the failures. However, only the 
lowest few percent of the population failure times will occur within 
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Fig. 3—The two lognormal curves (dashed lines) show possible relationships between 
accelerated-stress results and the basic failure rate model. 
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the service life. Therefore, accelerated testing does not accurately 
predict the failure rate distribution during the service life. For the 
lognormal examples shown in Fig. 3, only two percent of the compo- 
nents will have failed before the maximum failure rate is reached (in 
Fig. 1, at forty to four hundred years). 

Since neither direct observations nor accelerated life testing defines 
the time dependence of the long-term failure rate, how then do we 
justify using the constant-failure-rate exponential model? We have 
used the exponential model largely by default, since it is almost 
universally used to model long-term failure rates as, for example, in 
MIL-HDBK-217D. The model is a reasonable compromise between 
the decreasing failure rate of the Weibull model and the increasing 
failure rate of the lognormal curve. 

Field tracking does not define the failure rate during the very long 
term, but it can give failure rates at the end of infant mortality and 
the beginning of the long term. From those we can estimate the long- 
term failure rates. Accelerated life testing results are used in addition 
to tracking data. From these results we predict the total fraction of 
the population estimated to fail within the service life, which must be 
consistent with the failure rate estimates. If it is not, the estimates 
must be reevaluated. 


5.1 Sources of component failure rates 


If field tracking or accelerated testing data are available, the pro- 
cedure just described is used to estimate the long-term failure rate. 
Where we do not have any relevant data, we use MIL-HDBK-217D 
numbers. With few exceptions, our estimates for semiconductor com- 
ponents come from AT&T Bell Laboratories data and those for 
nonsemiconductors come from MIL-HDBK-217D. 

The infant mortality data that exist are solely from semiconductor 
components. These data form the basis for our estimates of the infant 
mortality parameters for all semiconductor components. 

We do not have good infant mortality data on nonsemiconductor 
components. We do, however, believe that infant mortality exists for 
these components. Because we do not have good data, we chose a value 
of a = 0.6 at the lower end of the range observed in equipment 
(between 0.6 and 0.9). Furthermore, we assume the infant mortality 
failure rate at 10,000 hours equals the long-term failure rate. These 
two assumptions, taken together, describe the existence of infant 
mortality, to a modest extent, in nonsemiconductor components. 


5.2 Effect of temperature 


Up to this point, we have only described the failure rate for a 
component as a function of time. In reality, a component’s operating 
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environment will also affect the failure rate. The component’s oper- 
ating temperature is one such environmental factor. It can have a 
strong effect on the failure rate. 

We base the component failure rate estimates on an assumed 40°C 
“typical” ambient temperature, since this is the temperature at which 
most of the available field tracking data are taken. Ambient temper- 
ature refers to the temperature in the immediate vicinity of the 
component. In most cases the failure rate estimates can be used 
directly. In cases where the temperature departs significantly from the 
typical value, its effect must be taken into account. 

The effect of temperature is usually modeled through its effect on 
the rate of aging of a component. At a higher temperature, failures 
will generally occur sooner and, therefore, at a greater rate. The 
difference in rates of component aging at two different temperatures 
is described by an acceleration factor. If T2 is a higher temperature 
than T;, then the Arrhenius relationship gives the following accelera- 
tion factor: 


A(T», T;) = oaln 3 


where kg is the Boltzmann constant and E, is the “activation energy”. 
T, and T> are in units of degrees Kelvin. This Arrhenius relationship 
is well understood for chemical reactions. However, its use in the 
current context is based purely on empirical evidence. Therefore, the 
constant E, does not really have physical meaning as an activation 
energy. Rather, it should be considered as an empirical curve-fitting 
constant. 

With the exponential model, a factor of A(T, T,) increase in the 
rate of aging leads to a factor of A(T, T)) increase in the constant 
failure rate, that is: 


AT, = A(T», T;)\r,- 


Under the Weibull model for infant mortality, the effect of tempera- 
ture on the failure rate is not as simple. A factor of A(T2, T;) increase 
in the rate of aging leads (after some careful algebra) to a factor of 
A(T», T;)'~? increase in the failure rate: 


Ar,(t) = A(T2, Ty)*"*A7,(t). 


It should be noted that the activation energy, and hence the acceler- 
ation factor, in the long term is not necessarily the same as during 
infant mortality. It depends on whether the expected cause of failures 
(failure mechanism) is the same. Table I lists the activation energies 
used. 
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Table |—Activation energies for selected components 


Component E, (eV) Reference 

Infant Mortality 

All components 0.4 4 
Long Term 

Discrete semiconductor 0.4 Unpublished work 

Bipolar integrated circuits 0.4 Unpublished work 

Metal-Oxide Semiconductor (MOS) inte- 0.5 Unpublished work 

grated circuits 

Ceramic capacitors 1.0 6 

Plastic capacitors—metallized and foil 0.12 6 

Film resistors—metal or carbon 0.08 6 

Carbon resistors 0.34 6 


5.3 Effect of electrical stress 


The level of electrical stress at which a component operates can 
also affect the failure rate. The higher the level of electrical stress, the 
more quickly we expect a component to fail. The effect of electrical 
stress, as for temperature, is modeled with an acceleration factor. This 
gives the difference in rate of aging at different values of applied 
electrical stress. The relationship we use to describe S, the acceleration 
factor, is 


S( po, D1) = emP2-Pr) 


where p2 and p, are stress levels. These stress levels are given as a 
percentage of the maximum specified level. The electrical parameters 
that constitute electrical stress are different for different types of 
components (see Table II). 

The stress constants (m) are based on information in MIL-HDBK- 
217D for long-term operation. There is no effect of electrical stress on 
integrated circuits, since the applied voltage is specified and assumed 
to be constant. Due to lack of information, no effects of electrical 
stress on failure rates during infant mortality are modeled. 

The level of applied electrical stress assumed for typical operation 
is 25 percent. If actual levels of applied electrical stress significantly 
differ from 25 percent, then the acceleration must be taken into 
account. As with temperature acceleration, the electrical stress accel- 
eration factor simply multiplies the constant long-term failure rate, as 
follows: 


Np, = S( po, P1)p,- 


Figure 4 shows how the combined effects of elevated operating 
temperature and high electrical stress can affect a component’s failure 
rate. Both temperature acceleration (Ayr) and acceleration due to 
electrical stress (S) can affect the long-term failure rate. Only temper- 
ature acceleration (Am) affects the failure rate during infant mortality. 
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Table 1l—Electrical stress dependence (5S) of 
selected components 


Electrical Stress Range of m 
Component Parameter Values 
Integrated circuits Not applicable — 
Resistors Power 0.006-0.024 
Capacitors Voltage 0.024—0.150 
Switches Current 0.013 


Note that these acceleration factors change the time at which the 
long-term failure rate is reached in the model. 


5.4 Effect of dynamic burn-in screening 


One widely used method to reduce the impact of infant mortality 
on equipment is the “screening” method. This refers to some activity 
performed on components or equipment to screen or “weed out” infant 
mortality failures before customer use. Components or equipment are 
stressed in some way and then tested. Any failures are removed or 
repaired. By inducing these failures to occur prior to use, customers 
should experience fewer equipment failures. Some commonly used 
screens are thermal cycling, high-voltage stress, or simply electrical 
operation (burn-in). 
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ae 4—Effect of temperature (Air and Am) and electrical stress (S) on failure rate 
model. 
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A variety of screens are believed to be effective in some cases. 
However, we have formally modeled the effects of only dynamic burn- 
in. This refers to the electrical operation of components or systems in 
a manner simulating eventual use. It consists of powering up and 
dynamically exercising the component or system for a period of time. 
This is distinguished from static burn-in in which power is applied 
but no dynamic exercise takes place. The burn-in may or may not 
occur at an elevated temperature. 

We model only dynamic burn-in because we do not have a suffi- 
ciently good understanding of the quantitative effects of other screens. 
Even the understanding of the effects of dynamic burn-in is poor. 
Mathematically, our burn-in model follows naturally from the basic 
failure rate model. However, there is little direct evidence to substan- 
tiate the model. 

One assumption provides the basis for the dynamic burn-in model. 
We assume that the failure rate of a component depends on the length 
of time of previous operation, wherever that operation occurred. As is 
clear from our model in Fig. 1, the more operating time a component 
has accumulated within the infant mortality period, the lower will be 
its failure rate. Therefore, a manufacturer can reduce the failure rate 
a customer will experience by operating components or equipment for 
a period of time before shipment. 

Calculating the effect of dynamic burn-in is a matter of calculating 
the effective operating time to which the burn-in is equivalent. The 
effective operating time, tes, corresponds to operation at the nominal 
40°C. Then, if t represents the amount of operating time after the 
burn-in, the failure rate at 40°C is 


A(t) = Arlt + tere). 


Operation after burn-in might occur at a temperature higher than 
40°C. In that case, the above failure rate is modified by the temperature 
acceleration factor already discussed. 

Figure 5 illustrates the modeled effect of dynamic burn-in. The 
straight dashed line shows the basic infant mortality failure rate. The 
solid curve gives the failure rate after a burn-in equivalent in time to 
ter. Note that this curve is a simple replotting of the dashed line but 
starting at age tec rather than at zero age. 

To calculate tes, we again make use of the temperature acceleration 
factor, Ar, introduced earlier. If burn-in occurs at a temperature 
above 40°C for x hours, then the effective burn-in time is ter = Arx. 
(Am is the acceleration factor, for the burn-in temperature, relative 
to 40°C.) 

Burn-in can be performed at any of several stages. Components can 
be burned in. Burn-in can also occur at the circuit pack or at the 
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Fig. 5—Effect of burn-in on early system failure rate. 


system level. Moreover, equipment can be burned in at more than one 
of these stages. Under the model, the effect of burn-in is cumulative. 
With more than one burn-in, the effective burn-in times at each stage 
are added to give the total effective burn-in time which is 

all stages 


(tets)total = p (tete)i. 


Caution should be exercised in using this additive assumption of 
burn-in treatments. In some instances the initial failure rate of a 
population of devices during a second burn-in was larger than the 
final failure rate during the first burn-in. This “setback” may have 
been due to the testing of the devices and insertion of the devices into 
circuit boards, which was done between the two burn-in treatments. 
This additional handling may damage some devices (e.g., by electro- 
static discharge) causing them to fail sooner than they would have 
otherwise. 


VI. FAILURE RATES FOR EQUIPMENT 


Up to this point, we have described the failure rate models applied 
to individual components. It is a simple matter to combine the com- 
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Table III—Environmental application factors 


Environment E 
Permanent structures, environmentally controlled 1.0 
Ground shelters, not temperature controlled 1.1 (Ref. 7) 
Manholes, poles 1.5 (assumed) 
Vehicular-mounted equipment 8.0 (Ref. 7) 


ponent failure rates to estimate the failure rates of equipment. Basi- 
cally, component failure rates are added together to give equipment 
failure rates. There is, however, one final modification we make to the 
summed failure rates. This modification accounts for a failure rate 
effect, which we do not understand well enough to apply at the 
component level. Rather, we apply it at the equipment level. 

The modification involves an environmental application factor, E. 
It reflects environmental factors other than temperature that affect 
the equipment failure rates. Values of E are listed in Table III. These 
cover the usual environments for AT&T telecommunications equip- 
ment. They are based on information in Ref. 7. 

With the inclusion of the equipment-level application factor, equa- 
tions giving equipment failure rates can be written. For infant mor- 
tality the equipment failure rate is 


all components 


Mott) = EY (Ami “(dt + (tere) 7%. 


The long-term failure rate is 


all components 


Atta = E y (Aur) SAL); 


i 


VII. PREDICTING THE NUMBER OF INFANT MORTALITY FAILURES IN 
A TIME INTERVAL 


Being able to predict the numbers of failures during the infant 
mortality period is important for a number of reasons. The estimates 
are useful for anticipating customer reaction. They can be used to 
understand warranty repair costs and to plan repair strategies. If any 
of these factors appear undesirable or unacceptable, the design, screen- 
ing, or requirements of the products can be reevaluated. If this is done 
at an early enough stage, changes can be made when the impact on 
cost is low. Predictions also provide standards against which compo- 
nent or equipment performance can be measured. Once production 
begins, the results of ongoing reliability testing can be compared to 
the standards to show where to concentrate ongoing efforts to make 
improvements. 

One calculation that is particularly useful is predicting the number 
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of failures out of some population during a stated time interval. The 
fraction of units failing between ¢, and tz is 


tg 
fraction fail = { A(t) dt, 
ty 
where X(t) is the unit’s time-dependent failure rate. This approach is 
strictly correct when all failures in a population are repaired. In 
addition, it is correct if they are replaced with units of the same age. 
However, it is also a good approximation even if failing units are 
replaced by units of different age or not at all, if the fraction failing is 
small (less than a few percent). 

The following example illustrates the strong effect of infant mor- 
tality. Typical failure rate parameters for integrated circuits are \, = 
10 FITs, a = 0.8, and X; = 16,000 FITs. (One FIT equals one failure 
per 10° component-hours.) In the long term, we would expect about 
10-° failures per component-hour, or 7.2 x 107° failures per compo- 
nent-month. In a system made up of 10,000 such components, we 
would expect 0.072 failures per system during one month. In the infant 
mortality period, the situation is drastically different. In the first 
month, we would expect 


720 
fraction fail = (1.6 x 107°) f iF dit. 
0 


or roughly 0.0003 failures per component. This gives three failures per 
system in the first month. The number of failures expected during the 
first month is roughly forty times higher than expected in one month 
in the long term. Clearly, such an effect, if not anticipated, could lead 
to nasty surprises. 


VIII. COMPARISON WITH MIL-HDBK-217D 


Many of the concepts embodied in our failure rate model are also 
used in MIL-HDBK-217D, the de facto industry standard for failure 
rate prediction. There are, however, some important differences be- 
tween the two models as well. As we already mentioned, the largest 
difference is in our formal inclusion of a model for infant mortality. 
This enhancement overcomes the major shortcoming of the MIL- 
HDBK-217D methodology, especially as applied to commercial prod- 
ucts. 

A second difference lies in the magnitudes of long-term failure rate 
estimates for integrated circuits. Our numbers are based on predives- 
titure Bell System experience. They are generally lower than those in 
MIL-HDBK-217D. Our data show that the MIL-HDBK-217D failure 
rates for integrated circuits are unrealistically high, especially those 
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for Large-Scale Integration (LSI) components. A final difference is in 
our dynamic burn-in model, which predicts an effect of burn-in on the 
infant mortality but not on long-term failure rates. MIL-HDBK-217D 
assumes that burn-in impacts long-term failure rates. 


IX. SUMMARY 


We have described the basic component reliability model used at 
AT&T Bell Laboratories to predict component and equipment relia- 
bility. It is a two-part model. A Weibull model describes infant 
mortality. An exponential model describes long-term behavior, beyond 
roughly the first year. The infant mortality part of the model is very 
important. It quantitatively describes the initially high, but rapidly 
decreasing, early-life failure rates. Recognizing such behavior is be- 
coming critical as increasingly complex electronic systems are being 
sold to a variety of customers. Widely recognized methods of predicting 
reliability, such as MIL-HDBK-217D, do not model infant mortality 
effects. 

Given the basic component reliability models and failure rate esti- 
mates, failure rates of equipment can be easily estimated. These 
equipment estimates, including infant mortality, are essential when 
planning for new, competitive product offerings and manufacturing 
them to have well-controlled reliability. 
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Waiting Time Convexity in the M/G/1 Queue 


By D. L. JAGERMAN* 
(Manuscript received June 8, 1984) 


Strong bounds are obtained on the complementary waiting time distribution 
for the M/G/1 queue using the a-convexity structural characteristic of the 
distribution. This notion is discussed and a sufficient condition is obtained. 


I. INTRODUCTION 


This paper investigates the a-convexity’ of the complementary 
waiting time distribution in the M/G/1 queue and shows how a 
sufficient condition for a-convexity is obtained, in terms of the service 
time distribution. This structure characteristic will permit strong 
bounds to be obtained on the complementary waiting time distribution. 
For convenience, the definition of a-convexity and some of its prop- 
erties are given below. 


Il. a-CONVEXITY 

A function f(x) is said to be a-convex on an interval I if e**f(x) is 
convex on J. Of course, ordinary convexity corresponds to a = 0. A 
sufficient condition for a-convexity is 


sae OE 
eo (efx) 20, x EL. (1) 


This is the same as 
f’(x) + 2af’(x) + a?f(x) = 0, xEl. (2) 


*AT&T Bell Laboratories. 
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A function may be a-convex without being convex; for example, 
consider f(x) = x°, which is a-convex for x = 0 and a = 0. For a = 1, 
however, x? is a-convex for -3—V3 < x < —3+¥V3 as seen by use of 
(2). 

The a-convexity of a function may permit stronger bounds than 
convexity to be obtained on the function or on integrals of the function. 
For example, let p(x) = 0 and let f(x) be convex on J; then Jensen’s 
inequality” states 


i f(x)p(x)dx = f(u) f p(x)dx, 


= i) xp (x)dx i | p(x)dx. (3) 


If f(x) is a-convex on J, then, since 


{ f(x)p(x)dx = f ef (x)e~*p(x)dx, (4) 


one has 


i, f(x)p(x)dx = e“*f(u) f e ““p(x)dx, 


= Jse~pemas / | empire (5) 


This result can be stronger than (3). 
An example is provided by 


) em 
ae J 1+ Pa (6) 


whose value is K = 0.5963. From (3), one has K’> 0.5. Since, for x = 
0, 1/(1 + x) is a-convex for all a, one may apply (5) to obtain 





a 





ae 
= 
at2’ 7) 
which, for a = (v5 — 1)/2, yields K > 0.5596. 
Let f(s) be the Laplace transform of a function f(x), that is, 
f(s) = i ef (x)dx, (8) 


and let the transform be absolutely convergent for s > 0. Then the 
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approximation sequence, f,(x) (n = 0, 1, 2, ---), introduced in the 
Laplace inversion theory, is given by’ 


fulx) = (oO) sh1FlON(s) | (9) 
n} pgs 
from which, in particular, 
Lng fl 
fo(x) = z. (2), 
fil) =- 7’ (2) | (10) 
x i 


The approximation sequence may be used to obtain bounds on f(x); 
thus, if f(x) is convex for x = 0, then 


f(x) <fnia(x) Sfrlx), x20, ned. (11) 


If f(x) is a-convex on x = 0, then the bound of (11) may be strength- 
ened. Let f(s — a) be absolutely convergent for s > 0; then it is the 
transform of a function g(x) for which 


f(x) = e g(x). (12) 
Application of (11) now provides the inequality 
f(x) S e “Snai(x) S e gn (x), x = 0, n= 0. (13) 


This is a much tighter inequality than (11), especially for the tail of 
f(x), and constitutes the main tool for bounding the M/G/1 waiting 
time distribution. 

The Bernstein theorem,*® which states that f(x) = 0 if f(s) is 
completely monotone and, conversely, may be used to translate con- 
dition (1) in terms of f(s). Thus let f”(x) be continuous on (0, «); then 
f(x) is a-convex on (0, ©) if and only if 


(s + a)*f(s) — (s + 2a)f(0+) — f’(0+) (14) 


is completely monotone in s on (0, ©) and is absolutely convergent for 
s>0. 

A function, f(x) > 0, is said to be log-convex if In f(x) is convex on 
some interval J. The condition for log-convexity is 


f" (x)F(x) — f’ (x)? = 0, xEl. (15) 


In particular, log-convexity implies convexity; hence e**f(x) is convex, 
and a log-convex function is a-convex for all a. The converse is also 
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true. This follows from (2) on observing that the discriminant of the 
quadratic in a is f’ (x)? — f”(x)f(x); hence a-convexity for all a implies 
(15) and the log-convexity of f(x). An interesting corollary of this is 
that the sum of log-convex functions is log-convex since, clearly, the 
sum of functions convex for the same a is again convex for this a. 
This theorem and eq. (14) permit ascertaining the log-convexity of 
f(x) from its Laplace transform. 


Ill. @-CONVEXITY IN M/G/1 


The starting point for this investigation of a-convexity in M/G/1 is 
the Pollaczek-Khintchine formula.* Let B(x) be the service time 
distribution and F(x) the complementary waiting time distribution; 
also let B(s), F'(s) be the corresponding Laplace-Stieltjes transforms. 
Then 


ps — Af{1 — B(s)] 

s— Al — B(s)] ’ 
in which ) is the arrival rate, » is the service rate, and p = \/p is the 
offered load. It is convenient to use the forward recurrence time 


distribution, 6(x), corresponding to B(x). Since the Laplace-Stieltjes 
transform, 6(s), of 6(x) is 


A 1 — B(s) 
io 


6(s) = ; (17) 
S 
one has 
se 1 — &(s) 
(8). = = 
(s) = p 1 — pils) 
F(s) = p eee (18) 
S$ 1— p@(s) 
in which F(s) is the corresponding Laplace transform. Clearly, 
A(s)~=, sm; (19) 
hence 
am 1- 
Bis) ~p-ME=P 5 see, (20) 
Thus, 
F(0+) =p, = F’(0+) = —pp(1 — p). (21) 


The information is now available to apply condition (14). That expres- 
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sion now takes the form 
N(s) 


1 pits)’ (22) 


2 
N(s) = up(1 — p) + — 
ap| « 
= {oa — p(s +A + 2a) + a A(s). 
The function 1/[1 — p6(s)] is the Laplace-Stieltjes transform of a 
monotone increasing function on (0, ©). If we write N(s) in the form 
2 
Ns) = a + (1 — p)B(s) 


a= pe +2) 1 _ By 


2 
-*S11 - B6)}. (23) 


we see that N(s) is a Laplace transform. For the function of (22) to 
be completely monotone, it is therefore sufficient that N(s) be the 
transform of a nonnegative function. If we let b(x) be the service-time 
density function, and r(x) the corresponding rate function, that is, 





_ d(x) 
r(x) = Ba) BG)’ (24) 
this condition may be written in the following two forms: 
of 25) soon B(x)] 
ee ot ee > BN = 
A(1 — p) 
a? = 
+ [1 — B(u)]du, 
1 — po 
2 _ 
ee ee Oe (25) 


u(1 — p) 1 — B(x) 


to assure the a-convexity of F(x). 
The case a = 0 of (25) yields the interesting result that convexity 
of F(x) is guaranteed by 


r(x) =X, x20, (26) 
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which also implies 
1— B(x) <e™, x20. (27) 


Application of (25) will now be made to the following class of 
distribution functions B(x): 


B(x) = i, (1 -— e ™)dG(u) (28) 


in which G(x) is a distribution function on (0, %). Condition (25) now 
becomes 





= res | 
e™lu—-rX—-—2at —]}dG(u) = 0. (29) 
0 1 u 


=p 


Now the integrand is nonnegative if 








? 4a? 
2u>rX+ 2at+ (A + 2a) =e : (30) 
— p 
which suggests the introduction of the quantity c defined by 
c = inf[x; G(x) > 0]. (31) 
Thus condition (30) need be satisfied only for u > c, and hence 
. 4a? 
2c >A+ 2at (A + 2a)" - = (82) 


assures (29) and the a-convexity of F(x). One implication of (32) is 
the following constraint on a: 


a <c(1 — p) — Ve(1 — p)(A — cp). (33) 
Application of (33) to the exponential case B(x) = e~” yields 
a < pu(1 — p), (34) 
which, of course, is consistent with the known result 
F(x) = pe~“t-e, (35) 


As another illustration, consider 
_ 1 —x 1 —2x 
B(x) =1 5° 3° (36) 
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for which c = 1, np = 4/3. One has 





3 1 3 
<1l--\-- —-)}. 
re eae a (1 ay (37) 
Since 
A i Pees 1 
B(s) == ania 
Oey 6 qo! (38) 
use of (18) yields 
ES ON 38s + 5 
EO ae (39) 
Sa Ae 2 = 5h 
If we let 
—3+rA+v1+ 
5 ei iii 
Pies Se 2 ea 
= 5 
Aa vit +14 3d 
= 5 
38V1+ 7-1-3) 
B= ————_-,, (40) 
2 
then we have 
F(x) = ees (Ae + Be*) (41) 
4AV1 + »? , 


This distribution is, in fact, log-convex; thus (37) is overly restrictive. 
This might have been expected since (25) is only a sufficient condition 
for a-convexity. 


IV. BOUNDS 


When values of a have been determined by use of the complete 
monotonicity of N(s) in (23), or by use of (25), then (13) may be 
applied to F(s) in (18). It is, of course, advantageous to use as large a 
value of a as possible consistent with the requirement that F(s — a) 
be absolutely convergent for s > 0. Applying (10) and (13) to (18) 
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Table I—Numerical inversion for different a 


F(x) 
Exact a= a= 
xX Results a=0 0.34549 0.69098 
1 0.1718 0.1815 0.1753 0.1727 
2 0.0834 0.0997 0.0883 0.0841 
3 0.0414 0.0590 0.0460 0.0416 
4 0.0207 0.0368 0.0245 0.0208 
5 0.0103 0.0239 0.0133 0.0104 
6 0.0052 0.0161 0.0073 0.0052 
7 0.0026 0.0112 0.0041 0.0026 
8 0.0013 0.0080 0.0023 0.0013 
9 0.0007 0.0058 0.0013 0.0007 


provides the following explicit bounds: 


ax 


e 





F(x) < 1- 


1 


ga a 
1-0 (7a 


As a numerical example, the approximation F(x) was calculated for 
F(s) of (39) with \ = 0.5. Table I shows the exact results obtained 
from (41). It also shows the inversion with a = 0 (that is, no a- 
enhancement used), the results with the value of a obtained from (37) 
(0.34549), and, finally, the calculation using the optimum choice, 
y = —a, i.e., a = 0.69098. 

As expected, the table shows improved accuracy as a is increased 
and, in particular, in the tracking of the tail behavior. This, of course, 
is the primary goal of a-enhancement. Observe that the approximate 
values are larger than the exact values, as implied by the a-convexity 
of F(x) and (13). 
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Accumulation of Jitter: A Stochastic Model 


By C. CHAMZAS* 
(Manuscript received March 29, 1984) 


A problem that has been considered extensively in the past is the accumu- 
lation of jitter in a chain of regenerative repeaters. For simplicity it is usually 
assumed that all repeaters in the chain are identical. However, this is not the 
case in real Systems, where considerable differences among the repeaters of a 
chain have been observed. These random variations are due mainly to manu- 
facturing tolerances, aging effects, temperature changes, etc. In this work, we 
examine the accumulation of systematic and random jitter along the chain, 
when the repeater transfer functions are subjected to random variations. We 
derive expressions for the expected value and variance of the output jitter 
spectrum in terms of average values of the repeater’s jitter transfer function. 
In addition we find the expected value and the variance of the RMS jitter. 
Finally we examine two special cases where the timing circuit employs either 
a phased-locked loop or a surface acoustic wave filter, and derive some 
asymptotic relations for the power spectral density and root-mean-square 
value of the accumulated systematic jitter. 


I. INTRODUCTION 


In a chain of self-timed regenerative repeaters for data transmission 
with Pulse Amplitude Modulation (PAM), each regenerator extracts 
the timing information (clock) directly from the received pulse train. 
Ideally the output of the timing circuit should be a sine wave with 
frequency equal to the baud rate of the data. In practice, however, 
imperfections of the circuits and noise in the transmission channel 
disturb the timing recovery operation so that the phase of the recovered 


* AT&T Bell Laboratories. 
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clock is randomly deviated from the desired input clock phase. Such 
deviations, called timing jitter, produce a position modulation of the 
regenerated signal; the timing jitter tends to accumulate along the 
chain and degrade the system’s performance. If the jitter introduced 
by the repeater’s timing recovery circuitry is the same for each repeater 
(for example, if it depends on the data pattern), it is called systematic 
jitter. Otherwise it is called random Jitter. 

The problem of jitter accumulation in a chain of regenerative 
repeaters has been considered extensively in the past.” For simplicity 
it is usually assumed that all repeaters in the chain are identical. 
However, this is not the case in real systems, where considerable 
differences have been observed among repeaters.® These random var- 
iations are mainly because of manufacturing tolerances, aging effects, 
temperature changes, etc. In this work we examine the accumulation 
of systematic and random jitter along the chain, when the repeater’s 
jitter transfer functions are subjected to random variations. Previous 
papers’* have treated the relationship between the jitter transfer 
function and the retiming circuit parameters. 

In this paper, we derive expressions for the expected value and the 
variance of the output jitter spectrum in terms of averages of jitter 
transfer functions of the ensemble of regenerators. In addition we find 
the expected value and the variance of the Root-Mean-Square (RMS) 
jitter. 

We also examine two special cases where the timing circuit employs 
either a Phased-Locked Loop (PLL) or a Surface Acoustic Wave 
(SAW) filter. Finally we consider the case of a retiming circuit exhib- 
iting random variations in its phase only and derive some asymptotic 
relations for the power spectral density and the RMS value of the 
accumulated systematic jitter. Of particular interest perhaps is relation 
(40), where we show that the RMS value of the accumulated systematic 
jitter for a nonpeaking case is approximately equal to V0.36*N/a, 
where N is the number of regenerators and a the phase slope of the 
average jitter transfer function. 


Il. THE THEORY 
2.1 The basic model 


A brief review of jitter accumulation is presented here for complete- 
ness. The model we use to study the jitter accumulation is similar to 
the linear model attributed to Chapman, as described by Byrne et al. 
Questions regarding the validity of this linear model will not be 
considered here. This problem will be addressed in a forthcoming 
paper. An oversimplified treatment is also given in Appendix B. 
According to this model, the timing jitter produced in a regenerative 
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Fig. 1—Chapman’s model for timing jitter accumulation in a chain of repeaters. 


section is the filtered sum of the jitter coming from the previous 
section plus an additional equivalent jitter e;(t) inserted at the local 
input. This is illustrated in Fig. 1, where W;(w) is the jitter transfer 
function of the ith repeater. The major differences between our model 
and the ones used previously are that we do not assume W;(w) to be 
the same for every repeater. 

The input jitter, e;(t), can be separated into its two components, the 
random part e;-(t) and the systematic part e;,(t), i.e., 


e;(t) = e,(t) + es(t). (1) 


The random component e;,(t) is different for each repeater and is 
usually caused by random sources, e.g., line noise, thermal noise, 
crosstalk, etc. The systematic component e;,(t) is the same for each 
repeater and is usually pattern dependent. 

Let &(w), &,(w), &;(w) be the two-sided power spectra of e;(t), e;-(t), 
and e;,(t), respectively. Assuming that e;,(t) and e;,(t) are statistically 
independent, then 


P(w) = &,(w) + (a). (2) 
If S(w) is the spectrum of the output jitter then we can write 
S(w) = ®,(w)T,(w) + ®,(w)T;(w), (3) 


where! 


T.(w) = | Wi(w) [7 + | Wi(w) We(w) |? 
+++) +|Wi(w) --- Wr(w)|? (4) 


is the total transfer function for random jitter, and 
T(w) = | Wi(w) + Wi (w) Weo(w) 
+ +++ + Wi(w)Wo(w) --- Www) |? (5) 


is the total transfer function for systematic jitter. 
If we assume that 


W;(w) = W(w) j=1,2,---,WN, (6) 


then we obtain the known relations? 
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1 — | W(w) |" 


T,(w) = | W(w) |? 1-|WwP (7) 
and 
»|1—- WX)? 
T;(w) = | W(w) | {=We) |" (8) 








As we will show, relations (7) and (8) provide us with a satisfactory 
approximation when the fine structure of S(w) is not important, as in 
the evaluation of the RMS jitter, 


ce = [ S(w)dw. (9) 


However, in some applications the detailed structure of S(w) is 
essential, and the random variations of W;(w) must be taken into 
consideration. Hence it is important to describe W;(w) with a proba- 
bilistic model instead of the deterministic approach used in the past. 

Another advantage of using this approach is that we can directly 
examine the impact of tolerances in manufacturing the repeater com- 
ponents on the accumulation of jitter. 


2.2 The stochastic model 


In this section we shall consider W;(w) to be a random variable, 
independent and identically distributed (i.i.d.). Thus S(w), the spec- 
trum of the jitter at the output of N regenerative repeaters, is also a 
random variable and it should be described with its expected value 
and variance. 

With the assumption that e;(t) is a stationary stochastic process we 
obtain from (2) and (3) 


E{S(w)} = ®,(w) E{T;,(w)} So ®,(w)E{T,(w)} (10a) 
Var{S(w)} = ©2(w)Var{T,(w)} + ©2(w)Var{T,(w)} 
+ 28,(w)®,(w)[E{T.(@)T-(w)} — E{T.(w)}E{T,(w)}], (10b) 


where E{X} and Var{X } are, respectively, the expected value and the 
variance of the random variable X. The spectra ®,(w) and ®,(w) of the 
input systematic and random jitter will be considered known, because 
we can measure them experimentally. Usually ®,(w) and ®,(w) are 
constant for the low frequencies where T,(w) and T,(w) are significant 
(white noise). Typical values for &, and ®, are between 1 and 100 deg?/ 
MHz. 

Thus for the evaluation of the E{S(w)} and Var{S(w)} in (10) 
we need to find E{T,(w)}, E{T,(w)}, Var{T.(w)}, Var{T,(w)} and 


E\T.(w)T,(w)}. The expressions for their analytical evaluation are 
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given below for a chain of N repeaters. Their derivations are given in 
Appendix A. For simplicity hereafter we shall omit the dependence of 
the various functions on the parameter w, except if it is not obviously 
implied. Since the derivation of Var{T,} and E{T,T,} is complicated, 
we simplified our analysis by making use of the Central Limit Theorem 
and assuming that the two random variables X and Y defined as 
X+jY =>", WW, --- W, are jointly normal. Hence relations (15) 
and (16) are not exact, because of the above approximation. 
For a chain of N repeaters 














= 1 — BN 
T, = E(T,| = B= (11) 
Ze 1— BN B 1- WY 1-B% 
= eit) = BF ore FG te 2) 
1— CN C 1-BN 1-CN 
2y _ Se ee 
pir = 0A So FB al oe) 
Var{T,) = E{T?} — E*{T;} (14) 


Var{T,} = 2[E7{X?} + E?{Y?} + 2E°{XY} — (E?{X} + E'{Y})?] (15) 
E{T,T,} — E{TJE{T,} = 2LE{(XJE(XT,} 
+ E{Y\E{YT,} — (E°{X} + E*{Y})E{T,}], (16) 


where 
W = E{W3} (17a) 
Z = E{W?} (17b) 
B= E{|W;|*} (17c) 
C = E{|W:\*} (17d) 
and 
D = E{|W;|?W3. (17e) 
With X and Y defined as 
x+jy=Q=3 W, We --- Wx (17f) 
we have 
E{X} = Re EQ}, ElY} = Im E{Q} (17g) 
E{X*} = [T, + Re E{Q*}]/2 (17h) 
E{Y*} = [T — Re E{Q*}]/2 (173) 
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E{XY} = Im E{Q"}/2 (17)) 














E{XT,} = Re E{QT;} (17k) 
E{YT,} = Im E{QT;}, (171) 
where 
Le 

EQ) = W-— (17m) 
» _ 71 -2N Z  fi1-w’ 1-2" 

Be = 24 +g 4 on) 
LD: D Law | pe 

Se ear ST +2 le Fl ed 


With the above formulas we can evaluate the E{S(w)} and Var{S(w)} 
in (10) in terms of the averages W(w), Z(w), B(w), C(w), and D(w) of 
the individual jitter transfer functions W;(w). These quantities can be 
estimated either experimentally from a sufficient number of samples 
of W;(w), or numerically, using a Monte Carlo technique, from an 
appropriate model of W;(w). The second approach will be used in 
our examples. If Wi(w) = W(w) i = 1, 2, ---, N then B = | W(w) |’, 
W = W(w), and after some algebra we can show that relation (12) is 
equivalent to relation (8). . 

Another quantity, beyond S(w), which is important in jitter accu- 
mulation, is o,, the RMS value of the output jitter. We shall also 
evaluate its expected value and its variance. 

The expected value is obtained easily from (9) and (10a) as 


oO 


Fe eee { E{S(o)}de. (18) 
27 J-x 

The evaluation of the variance of o2 is more difficult. This is done 
in Appendix A. However, the result is complicated because it involves 
the evaluation of R(u, v) = E{W,(u)W;(v)}, the autocorrelation of 
W;(w). To avoid the additional computations needed for the calculation 
of R(u, v), we obtain an upper and lower bound for Var{c2} by assuming 
that W;(w;) and W;(ws) are, respectively, highly correlated or uncor- 
related. Then we can derive (see Appendix A) that 


1 eo 1/2 
_ | { Var{S(w)}dw 
on —0 


where Var{S(w)} is obtained from (10b). The average value of the two 
bounds appears to be a good estimator of the true variance of the 
RMS. With the relations given above we can now evaluate E{S(w)}, 
Var{S(w)}, E{o2}, Var{o2}. 


ioe) 


< Var{a3} ei Var{S(w)}dw, (19) 
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In deriving E{S(w)} in (10a) and Var{S(w)} in (10b), we have made 
the additional assumption that (w), the spectrum of the additive 
jitter e;(t), is the same for all the regenerators. This appears to 
contradict our assumption of different W;(w), since it is well known 
that ®(w) is highly correlated with the transfer function of the re- 
generator. This problem can be resolved if we define e,,(t), the system- 
atic component of e;(t), as the part that is identical, within a constant, 
for all the repeaters, i.e., e;,(t) = a;e,(t), where a; depends on W;(w). 
The random component e;-(t), the remaining part of e;(t), has a 
spectrum ®;(w) and is statistically independent for each repeater. Then 
in relation (5) we have to replace W;(w) by a;W;(w), while in relation 
(4) we must replace | W;(w) |? by ®;(w) | W;(w) |”. Using these substi- 
tutions, all the remaining expressions are still valid. The evaluation 
of a; and ®;(w) in terms of the transfer function W;(w) is a difficult 
problem directly related to the problem of expressing the statistical 
properties of e;(t) in terms of the regenerator’s transfer function. This 
question has been considered in Refs. 4 through 6. 


I. APPLICATIONS 


In the following three sections we shall present some numerical 
applications of the above theory. In Section 3.1 we shall examine 
chains employing phase-locked loops in their timing circuits.*® In 
Section 3.2 the timing circuit employs a SAW filter.”® Finally, in 
Section 3.3 we shall examine a special case, where W;(w) is assumed 
to exhibit a random variation only in its phase while its amplitude is 
identical for all repeaters. This case is of interest when we have a 
maximally flat filter and we want to examine the spectrum of the 
accumulated jitter for low frequencies. 

All the numerical results are normalized by assuming ®,(w) = 
&,(w) = 1 deg?/MHz. 


3.1 Timing circuits with PLL 


Phase-locked loops have been used extensively in the timing circuits 
of regenerative repeaters and their jitter performance has been ex- 
amined thoroughly.*** If the timing extractor employs a PLL, its jitter 
transfer function W(w) is equal to the phase transfer function of the 
PLL, under the assumption that the output of the phase detector is 
small. This is equivalent to the assumption of having small alignment 
jitter. (Alignment jitter refers to the deviations in alignment between 
the clock embedded in the incoming data stream of the regenerator 
and the timing clock derived from the data stream by the timing circuit 
of the regenerator.) The above requirement is usually satisfied, vali- 
dating the linear model of the PLL. The phase transfer function of a 
PLL used in the literature is typically described with a second-order 
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model. If we want to consider parasitic elements in the PLL, then, as 
M. W. Hall suggested,’ a fifth-order model is more realistic, especially 
for high frequencies. However, the spectrum of the accumulated jitter 
is significant only in the low frequencies. Hence the second-order 
model provides an adequate description of the PLL in this frequency 
range, because it coincides with the fifth-order model in this range of 
frequencies. 

Thus the jitter transfer function W;(w) of the ith repeater can be 
modeled as 


2 
26; WniS = Whi 


Wi(s/j) = sea On 
(s/7) 82 + Wiwnis + wi 


$= ju, (20) 
where ¢; is the damping coefficient of the PLL and w,; its natural 
frequency. 

We shall assume that the average values of ¢; and w,; are ¢; = 6 and 
@nj = 4.5 kHz. Also, the bandwidth of the PLL is assumed to vary 
from 20 kHz up to 80 kHz. To simulate the above conditions we 
assumed that ¢ is uniformly distributed from 4 to 8 and w,; is also 
uniformly distributed from 2 kHz to 7 kHz. Using the above numbers, 
we calculated W(w), Z(w), B(w), C(w), D(w) (see 17a, b, c, d, e) by 
numerically averaging 1000 samples of W;(w). In Fig. 2 we show the 
amplitude and phase of W(w) as well as two samples of W;(w) with 
¢; = 6 and w,; = 2 kHz and 7 kHz. The shaded area shows the 
permissible range of W;(w). Bandwidth varies between 16 kHz and 
110 kHz, while phase slope varies between —0.6 deg/kHz and —3.2 
deg/kHz. In Fig. 3a we plot E{T,(w)}, and E{T,(w)} + Var{T,(w)} for a 
chain of 50 repeaters. Also in the same graph we show the 7;,(w) 
evaluated by using relation (8), i.e., assuming that all repeaters have 
identical PLLs with W(w) = E{W;(w)}. In Fig. 3b we show the same 
results for 200 repeaters. Figure 4 shows the same curves for the 
random jitter component. 

To check our theoretical results, we did a complete numerical 
Monte-Carlo simulation of 200 chains and computed the various 
parameters, including E{o2} and Var{c2}. The resulting graphs were 
indistinguishable from those shown in Figs. 3 and 4; the numerical 
values obtained for E{o2} and Var{c2} are shown in Table I. All the 
numbers agree with the results of our theoretical analysis. Also, the 
expected value of the RMS of the accumulated jitter, systematic and 
random, is shown in Fig. 5 with the +3 Var{o2} curves (99.7 percent 
confidence). Also in the same graph we plot the RMS values for the 
accumulated jitter for a chain having identical repeaters with jitter 
transfer function equal to'the average value E{ W;(w)}. 

We can draw the following conclusions from the previous example: 
(1) The model of identical repeaters having as jitter transfer function 
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Fig. 3—Total transfer function T,(w) for systematic jitter for (a) 50 repeaters and (b) 
200 repeaters with PLL using the stochastic model and the model with identical 
repeaters. 


the average jitter transfer function underestimates slightly the average 
RMS value and the power spectrum density. (2) The variance of the 
RMS jitter was only 5 percent (see Table I and Fig. 5) even if the 
bandwidth of the PLL jitter transfer function varied from 16 kHz up 
to 110 kHz. This implies that the model of identical repeaters provides 
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Fig. 4—Total transfer function T,(w) for random jitter for (a) 50 repeaters and (b) 
200 repeaters with PLL using the stochastic model and the model with identical 
repeaters. 


us with a reliable estimation of the RMS jitter. (3) The total transfer 
function T,(w) has a smaller variance than T',(w). 


3.2 Timing circuits with SAW filters 


The introduction of optical-fiber transmission systems has made 
possible data transmission rates of several hundred megabits per 
second. This introduced significant changes in the construction of the 
timing extracting circuits in the regenerative repeaters. The popular 
PLLs had to be replaced, because their implementation above 100 
Mb/s has been difficult, especially in integrated circuit form. Cur- 
rently, SAW filters’ have emerged as their replacements. This actually 
represents a return to passive filtering after many years of using the 
PLL. 

In this section we analyze a system where the tuned filter in the 
timing circuit is a passive SAW filter. In the present analysis we will 
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Table |—Theoretical and numerical values for E{o2} and Var{os} for chains with 50, 100, 200, and 300 repeaters 























E{o2} deg? * Var{os} 
Using Random Using Identical Numerical Result Numerical Result 
Model Repeaters (200 chains) Lower Bound Upper Bound (200 chains) 
N Syst. Ran. Syst. Ran. Syst. Ran. Syst. Ran. Syst. Ran. Syst. Ran. 
50 PLL 18.57 1.44 18.20 1.36 18.61 1.45 0.18 0.01 1.85 0.18 1.06 0.15 
SAW 28.94 8.45 28.82 8.00 0.08 0.05 0.93 0.90 
100 PLL 51.49 2.70 50.57 2.90 51.57 2.72 0.56 0.01 4.50 0.26 2.27 0.18 
SAW 58.17 13.94 57.98 13.01 0.23 0.08 1.53 1.47 
200 PLL 188.85 7.73 184.55 7.05 189.16 7.77 2.68 0.05 17.82 0.60 8.61 0.36 
SAW 116.67 22.85 116.38 20.97 1.02 0.11 3.27 2.18 
300 PLL 540.64 22.41 522.45 20.12 550.53 22.68 10.12 0.21 62.89 1.80 32.93 1.25 
SAW 175.19 30.43 174.83 27.58 2.02 0.15 4.87 2.80 


* The input jitter is ®,(w) = &,(w) = 1 deg?/MHz. 
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Fig. 5—RMS values of systematic and random jitter for repeaters with PLLs. Curves 
show: (1) stochastic model; (2) conventional model with identical repeaters; and (3) 
99.7-percent confidence interval. 


not consider the effects introduced by the prefilter and nonlinear 
device, which often precede the tuned circuit. It can be shown that 
this simplification is not restrictive.” A simple proof is also presented 
in Appendix B. 

Let H(w) = A(w)e’* be the transfer function of the SAW filter. 
Then, as has been shown,” the jitter transfer function of the regener- 
ator is approximately given by 


he wo el?) + Hw + wo)e leo 
2A (ao) ; 


where wy is the baud rate of the data. Relation (21) is valid under the 
assumption that the accumulated jitter does not have large compo- 
nents at high frequencies. The exact conditions for the validity of (21) 
will be studied in a forthcoming work. A simple derivation of (21) is 
given also in Appendix B. 
Let us now define the normalized low-pass equivalent of H(w) as 
H(w 4 wy )e 10) 


Hy,(w) = = ior for w > —a». (22) 


W (w) for |wl<q@, (21) 
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Thus W(w) is now given by 


_ [Hi (w) + Hi(-«)| 


W(w) 5 


(23) 
where « denotes complex conjugate. 

A model for transversal SAW filters related to its design electrical 
characteristics, i.e., number of fingers, distance of transducers, termi- 
nating impedances, etc, is given in Ref. 8. Since in the present 
simulation we are interested only in its transfer function, we will use 
a simpler representation, and model H;,(w), the low-pass equivalent of 
the SAW filter, as a two-pole filter. Its transfer function is 
$1 S82 


H,(s) = 


~ (s — 8:)(s — 82) 


e*, (24) 


where 
s; = —a; + jb; = —(1 + ea + Jb; as 


are its two poles and c is a linear phase slope used to correct the phase 
of the model, since SAW filters are not minimum phase filters. For 
jitter studies our model provides an adequate description for SAW 
filters. 

Denoting by BW the bandwidth of the SAW filter, w,. its center 
frequency, and wy the baud rate of the data, the various SAW filter 
parameters can be defined in terms of a, b;, and the asymmetry factor 
« as follows: 


We — Wo = (Bb; + b.)/2 mistuning (25) 
b = (b; - be) /2 a = (a; + a2)/2 
a = 2(@ — wo)/BW detuning factor 


wr =a’ + 0 natural frequency 
6 = a/op damping factor 
Q = w./BW filter’s quality factor. (26) 


The natural frequency w, can be determined approximately from 
the BW (for « « 1) by 


1) Gi ere ee 
= ( 2 7-3 +a eye 


For the simulation we used the following numerical values for the 
range of bandwidth BW, static offset, damping factor ¢, and asymmetry 
factor e: 
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160 kHz < = < 240 kHz 


lf. — fol < 50 kHz 
0.60 < ¢< 0.80 
-0.1<¢«<0.1 
c = —0.2 deg/kHz. (28) 


We assumed the above parameters to be uniformly distributed 
between their upper and lower limits. If the baud rate of the data is 
300 MHz, then the above values correspond to filters having Q’s 
between 625 and 940. To illustrate the relation of the above parameters 
to H,(w) and W(w), we plot in Fig. 6 the passband of the SAW filter 
and the corresponding jitter transfer functions for ¢ = 0.65 and various 
combinations of « and a. Notice that the frequency scale has been 
normalized to BW/2, the bandwidth of H;(w). 

In Fig. 6a we plot the passband of the SAW filters for ¢ = 0.65 and 
« = 0, 0.1. In Fig. 6b we plot the corresponding jitter transfer function 


AMPLITUDE IN DECIBELS 





1.0 -0.5 0 0.5 1.0 -1.0 -0.5 0 0.5 1.0 
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Fig. 6—(a) Underdamped SAW filter with ¢ = 0.65 and asymmetry factor « = 0 and 


0.1. Jitter transfer functions obtained with various detuning factors a = 2(a. — w.)/BW 
for (b) aw = 0 (0.105 dB jitter peaking), (c) a = —0.50, and (d) a = 0.50. 
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W(w) when f. = fo(a = 0). For this case the asymmetry is almost 
canceled. This is expected because from relation (23) we can see that 
asymmetries in H;(w) that are odd with respect to f, will be canceled 
when we form W(w). In Fig. 6c we plot the corresponding jitter transfer 
function when the detuning factor is a = —0.5, ie., the baud rate is 
100 kHz less than the center frequency of the filter. Figure 6d shows 
the jitter transfer function for a = 0.5. In Fig. 7 we plot E{W;(w)} 
(1000 samples) as well as two extreme samples of W;(w). Most of the 
W;(w) will be between W,(w) and W2(w). The shaded area shows the 
permissible range of W;(w). In Fig. 8 we plot the average total transfer 
function for systematic jitter, E{T,(w)}, and E{T,(w)} + Var{T,(w)} for 
50 and 200 repeaters. The variance of T,(w) in this example is much 
less than the variance of the example with PLL. This is due to the 
narrow distribution of the SAW filter phase slopes (—0.5 deg/kHz to 
—0.8 deg/kHz). Finally, in Fig. 9 we plot the expected value of the 
RMS of the accumulated jitter, systematic and random, with the 
+3 Var{od} curves (99.7 percent confidence). Some numerical values 
-for the accumulated RMS jitter and its variance are shown in Table I 
for N = 50, 100, 200, and 300. To find the true RMS value we have to 
multiply the numbers shown in Table I with the #,(w) and ®,(w) 
measured in deg”/MHz. For example, if N = 200 and a PLL is used, 
then assuming ®, = 20 deg?/MHz and 9%, = 5 deg”/MHz, we obtain 


o, = [20-189 + 5.7.73] = 61.8 degrees. (29) 
Using the same numbers for the SAW filters we obtain 
oy = [20-117 + 5-22.9]/? = 49.5 degrees. (30) 


Thus the simulated SAW filters accumulate less jitter, even with a 
bandwidth larger than the bandwidth of PLLs. This is due to the fact 
that PLLs have an inherent jitter peaking and because their phase is 
smaller than the SAW filters [see eq. (40)]. To obtain the correspond- 
ing peak-to-peak values we usually multiply the RMS value with a 
peak-to-peak/RMS factor. Typical values for this factor are between 
8 and 15. 

From the above example it becomes clear that, using our random 
model, we can tolerate larger manufacturing variances, because we 
can accept retiming circuits exhibiting substantial jitter peaking. For 
example, W,(w) (0.4 dB jitter peaking) in Fig. 7 can be accepted if the 
expected average jitter transfer function W(w) of the manufactured 
SAW filters possess a moderate jitter peaking (i.e., less than 0.1 dB). 


3.3 Timing circuits with random phase 


In this section we will consider the dependence of the jitter power 
spectrum on phase variations. This case is of interest when we want 
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Fig. 7—(a) Amplitude and (b) phase of E{W;(w)}, Wi(w), We(w) for the SAW model. 
W,(w) is obtained with BW/2 = 240 kHz, a = 0, £ = 0.60, and « = 0.1 (0.4 dB jitter 
peaking). W.(w) is obtained with BW/2 = 160 kHz, a = 0.25, § = 0.60, and « = 0.1. Area 
with lines indicates the permissible range of W,(w). 
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Fig. 8—Total systematic jitter transfer function T,(w) and its variance for 50 and 
200 repeaters. 
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Fig. 9—RMS values of systematic and random jitter for repeaters with SAW filters. 


to estimate the jitter spectral density for low frequencies.? We shall 
assume that 


W;(w) = A(w)eieo)*nlo)t (31) 


where 7;(w) are random variables with zero mean and independent 
and identically distributed (i.id.). Let ®,(X) be the characteristic 
density function of 7;, i.e., 


@,(X, w) = Efe *} = { f,(X, w)e"*dn. (32) 
Then [see relations (17)] all the needed statistics for W;(w) can be 
obtained analytically if &,(X, w) is known, ie., 
W(w) = E{W;(w)} = A(w)e* &, (1, w) 
Z(w) = E{Wi(w)} = A? (w)e”* ®, (2, w) 
D(w) = E{| Wi(w) |? Wi(w)} = A*(w)e?, (1, w). (33) 


If n:(w) is assumed to have zero mean and to be uniformly distributed 
between —W(w) and V(w), then 
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&,(X, o) = a 


If ;(w) is assumed to be zero mean and Gaussian distributed with 
variance o(w), then 


(34) 


®,(X, w) = ow? (35) 


Using (84), (35), (83), and (10) we can obtain E{T,(w)} and Var{T,(w)}. 
In Fig. 10 we plot E{T,(w)} for a transversal filter with maximally flat 
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Fig. 10—Systematic jitter power density versus phase for identical repeaters, repeat- 


ers with uniform phase distribution, and repeaters with Gaussian phase distribution for 
(a) N = 100 repeaters and (b) N = 180 repeaters. 
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transfer function (¢ = 2/2) as a function of its phase ¢(w) for 100 
and 180 regenerators for the cases n; = 0, n; uniformly distributed with 
W;(w) = 8.1° and n; gaussian distributed with o;(w) = 8.1°/3. 

In Fig. 11 we plot E{T,} and its variance for the case of uniform 
distribution. 

Finally, we would like to note the following approximate relations. 

1. Large variance approximation—lIf the variance of the phase is 
large for w = wo, then W(wo) « 1 because ®,(1, wo) is small [see 
relations (33) and (34)], and from relation (12) we obtain 


50 


we } - Var {75} 


-E{Ts} + Var {Ts} 
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Fig. 11—Mean and variance of systematic jitter power density versus phase for 
repeaters with uniform phase distribution for (a) N = 100 repeaters and (b) N = 180 
repeaters. 
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lim E{T;(wo)} = E{T,(wo)}. (36) 


a(wo)->% 
This is expected because the randomness of the phase removes the 
coherent accumulation of the systematic jitter for w = wo. 
2. Low frequencies and long chain approximation—For large N the 


jitter energy is concentrated near w = 0. For this region, we can assume 
that A(w) = 1 and approximate relation (12) with 


2[1 - W’]W — N(1 — W?) 
(1— wy w = 0. (37) 
3. High frequencies and long chain—For high frequencies where 
A(w) « 1 and large N) we obtain that E{T,} ~ E{T,} = E{| W;|7}. 
4. Small-phase variance approximation—If we assume that the 


variance o(w) of the phase ¢(w) is small, after some algebra we can 
obtain the following relation: 


E\T,} ~ Re \ 








N 
EXT;} = An t+ a” (w) Y An, (38) 
n=1 
where 
tic A” (wei) 2 
= De ee ee 
A, (w) — | A(w) | f= A(w)e??) 
and 
a7(w) < x. 


For low frequencies, w ~ 0, we can assume A(w) ~ 1 and the term 
A, (w) in relation (38) can be approximated with 


. now) | ? 
sin — 


Alo) = $b (39) 
sin “2 


5. Linear phase, low frequencies, small variance approximation— 
Let us assume ¢(w) = — aw for w = 0..Then o7(w) = S’w, where S is 
the variance of a, and relation (39) becomes 


. naw |? 
sin —— 





A,(w) = = 
sin 9 


For large N and assuming no jitter peaking and small-phase variance 
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we can approximate the variance o2 of the systematic jitter at the 
output of N regenerators with 


2 3 
E{o2} = ,(0) p= N+ - (“2*) ni deg”, (40) 


where 

®,(0) is the spectral density in deg’?/MHz 

a is the average phase slope in deg/kHz 

S? is the variance of a 

N _ isthe number of regenerators. 

Relation (40) shows that for long chains of regenerators exhibiting 
no jitter peaking the RMS value of the systematic accumulated jitter 
is determined mainly by the dc phase slope of the jitter transfer 
function and not by its shape. This is illustrated in Fig. 12, where we 
evaluate the RMS value of the systematic accumulated jitter for chain, 
consisting of identical regenerators having the jitter transfer func- 
tion the W(w) shown in Fig. 7. A curve using numerical integration 
is compared versus the curve predicted by the simple formula 
oo = 0.6V.N/a (a = 0.63 deg/kHz). For a first-order filter with band- 


~0.6/ N/a 


@ = 0.63 DEG/kHz 


DEGREES 
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Fig. 12—RMS systematic jitter accumulation using numerical integration and rela- 
tion (40). 
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width B, we have B = 0.36/a and the above relation becomes o, = 
VBN. This is the widely used relation (69) given in Ref. 1, which is 
valid only for first-order filters. 


IV. CONCLUSION 


We have presented a generalized model for the accumulation of 
jitter. This model differs from those used previously in that we do not 
assume all the repeaters have the same transfer functions. We have 
derived analytical expressions for the variance and the mean of the 
accumulated jitter in terms of the repeater’s jitter transfer function 
W;(w). 

We have also presented some numerical results for PLL and SAW 
filters. From our numerical simulations we found that for long chains 
the variance of the RMS jitter is about 5 percent. This implies that 
we can still reliably estimate RMS jitter by assuming that the jitter 
transfer function of all the repeaters is equal to the average jitter 
transfer function. Another result obtained from our modeling is that 
retiming circuits exhibiting large jitter peaking are acceptable if their 
average jitter transfer function does not have jitter peaking. 
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APPENDIX A 
Stochastic Evaluation of Jitter Parameters 

For simplicity the dependence on w is omitted in most of the 
following formulas. 
A.1 Evaluation of E{T,(w)} 

The transfer function T,(w) for the random jitter is given by [see 
(4)] 

T,=|Wil? +[|WiWol? + --- +|WiWe--- Wal? = (41) 
With the assumption that the W; are i.id., we obtain the expected 
value of T, as 
1 — BY 


E{T} =B+B?+Be+...+BY=B oan ye 





(42) 
where 

B= E{| W;|?} 
and this is the desired relation (11). 


A.2 Evaluation of E{T,(w)} 


The transfer function T,(w) for the systematic jitter is given [see 
(5)] by 











2 
W. --> Wh (43) 
The expected value of T, may be written 
N N 
k=1 m=1 
N 
-2{5 | W, Ws a iat W, |? 
k=1 
N k-l 
+2Re ¥ ¥ |WiW.e --- Wal? Wer Wise -°> Wp 
k=1 m=1 
Defining 
W=E{W}, B= E{|Wil7} (44) 
and using that W; are i.i.d. we obtain 
N k-1 
E{T,} = 5 Be +2Re > Y B™wr™ 
k=1 m=1 
— BN ,-1 1 — (B/W)F" 
= B r+ 2Re D BW am 
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or 





1 — BN _ B 1-—- WY 1-B% 
min} = BA sore _ 


which is the desired relation (12). 
For large w where | W;(w)| < 1 we can obtain from (42) and (45) 
that 


E{T;} = E{T,} = B(o). (46) 


A.3 Evaluation of Var{T,} 


To evaluate the variance of 7T;(w), the transfer function of the 
random jitter, we need only to find E{T?}, since 


Var{T;} = (E{T?} — E°{T;}) 
and E{T,} has been evaluated in (42): 


. N N 
E{T;} = E 12 [WiW.--» Wil? 2» |WiWe--- Waltt. (47) 
=1 m=1 
/ 
Defining B(w) as in (44) and 


C(w) = E{| Wi(w) | 4}, (48) 
we obtain 
N 


E{T?} -2{5 | W; We ses Ww, |* 
k=1 


N N-1 
+22 »y | Wi WoW * | Winer m+2 °°° walt 
=1 m=1 
N N N-1 
=yYcr+2¥ ¥ cmpem 
k=l k=1 m=1 


or 














1-c% Cc 1-BY 1-C% 
2 oe us 
pyr) = CAE 4g SEB Ta sy 


and this is relation (13). 


A.4 Evaluation of Var{T,(w)} 


A direct evaluation of the Var{T,} is possible but the resulting 
formula is lengthy. To simplify our analysis we will make an additional 
assumption. 

Let us define [see (17f)] 


‘N 
Q=) W.W.---Wr,=X+ iY (50) 
k=1 
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and 
X = E{X}, Y = E{Y}. 


Then, if N is large, since W; are independent we can assume that X 
and Y are jointly normal. This assumption will permit us to avoid 
calculations of fourth-order statistics. 

Since X, Y are jointly normal, the following relations are valid (see 


Ref. 10, p. 374): 
E{(X"} = 3E°2{X%} — 2X* 
E{Y*} = 3E7{Y?} — 2¥4 
E{X?Y?} = QE{XY} — 2X?Y? + E{X*}E{Y}. (51) 
Then 
Var{T,) = E(T?} - E{T,} 
= E{X? + Y’)} — EY{T,} 
= E{X*} + E{Y*} + 2E{X°?Y4} — E°{7)} 
= 3E7{X?} + 8E7{Y?} — 2x4 — 2Y? 
+ AE{XY} — 4X? Y? + QE{XAEY?} 
— E4X?} — E4Y4) — QE{X1E{Y?} 
or 
Var{T,} = 2{E°{X?} + E4Y?} + QE{XY}} — 2(K2 + ¥%)2, (52) 


which is relation (15). 
To evaluate (52) we need X, Y, E{X7}, E{Y*}, E{XY}. We evaluate 


these terms below. 
Since Q= X+ JY we have 


E{Q?} = E{X?} — E{Y?} + QjE{xy} 
E\T.} = E{|Q|7} = E{X*} + E{Y*}. 
Thus 
X = Re E{Q} 
Y = Im E{Q} 
E\{X*} = {E{T.} + Re E{Q?}}/2 
E\Y*} = {E{T.} — Re E{Q’}}/2 
E{XY} = Im E{Q*}/2. (53) 
Hence, we need evaluate only E{Q}, E{Q7}, since E{T,} has already 
been calculated in (45). From (50) and (44) we obtain 
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1— WN 
1-W- 


Also, following a similar method with the evaluation of E{T,}, we 


can find that 
1-2" Zz f1-WN 1-25 
ee, ee poe 7 ae. (55) 





E{Q} = W (54) 


EQ’) = Z 











where 
Z(w) = E{W?(w)} 
and this relation completes the evaluation of Var{T,}. 


A.5 Evaluation of E{T,T,} 


A direct evaluation of E{T,T,} is possible but the derivation and the 
resulted formula are lengthy. To facilitate our analysis, let us assume 
that Q = X + /Y and T, are jointly normal. Then 


E{T,T,} = E{X?T, + Y?T,}. (56) 
But since X, Y and T, are assumed jointly normal we have 
E{X?T,} = 2XE{XT,} — 2X°E{T,} + E{X°}E{T, 
E{Y°T,} = 2YE{YT,} — 2Y°E{T,} + E{Y7}E{T,} (57) 
and 
E\T.T,} — E{TJET;} 
= 2{XE{XT,} + YE{YT,} — E{T,}(X? + Y’)}. (58) 


All the above terms have been evaluated in (53) except for E{XT,} 
and E{YT,}. Since Q = X + jY and T, is real, it is enough to find 
E{QT,}. From relations (50) and (4) we obtain 


E\QT,} -F1 3 Wi W/W Walt 
k 1 


mm 





N K N m-1 
= b} b} D"™ wr + > > Di wk 
k=1 m=1 m=2 k=1 
or 
1 — DN D 1—-W* 1-—D)% 
elo ama Ta ete _ "|. (59) 
and then 


E{XT,} = Re E{QT;,} 
E{YT,} = Im E{QT;}. (60) 
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A.6 Evaluation of Var{c2} 
We have defined 


1 


cs on [ S(w)dw. (61) 


Hence 


Var*{o2} = Efot} — E{o2} 


= 2: [ i [E{S(u)S(v)} — E{S(u)}E{S(v)}]dudv 


= { { I(u, v)dudbv. (62) 
An” J-« 


To facilitate our calculations we will assume that the random term is 
small compared to the systematic term and we shall approximate the 
spectrum S(w) as 


S(w) = &,(w)T,(w). (63) 
A direct evaluation of Var{o?} is possible but requires knowledge of 
R,(u, v) = E{W;(u) W;(v)} 
R2(u, v) = E{W;(u) W; (v)}. (64) 


A simpler approach is to evaluate an upper and lower bound for 
Var{o3}. We can obtain an upper bound if we assume that S(u) and 
S(v) are highly correlated, i.e., 


E{(S(u) — S(u))(S(v) — S(v))} = Var{S(u)}Var{S(v)}, (65) 


and a lower bound if we assume that S(u) and S(v) are uncorrelated, 
1.€., 


E{(S(u) — S(u))(S(v) — S(v))} = Var{S(u)}5(u — v). (66) 
Then 


1/2 


ive] 


< Var(a3} < = i Var{S(w)}dw. (67) 





= | { Var{S (w)}dw 


Var{S(w)} has already been evaluated in (10b). It is expected that for 
long chains (large N) the true value will be closer to the lower bound, 
while for short chains (small N), it will be closer to the upper bound. 
From our numerical simulations we found the average of the two 
bounds to be a good estimator for Var{o?}. 
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APPENDIX B 
On the Jitter Transfer Function of a Tuned Circuit 


In most applications the jitter transfer function of a timing recovery 
circuit is approximated with the phase transfer function of the timing 
passive bandpass filter. In this appendix we derive in a simple way the 
phase transfer function of an arbitrary filter. The limits of the appli- 
cability of the derived formula are discussed. Finally, the results are 
extended for the case when a prefilter followed by a squarer is used. 


B.1 Phase transfer function 


In this part the phase transfer function of a narrow passive bandpass 
filter is shown to be 


H(w — woe’ + H(w + woe 


Wears 2A(wo) 


, for |w| <wo, (68) 


where H(w) = A(w)e’* is the transfer function of the bandpass filter 
and wo is the baud rate of the received data. In Fig. 13 we illustrate 
the above relation. 


H(w) = A(w)e/ Ol) 
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Fig. 13—Construction of the jitter transfer function W(w) from H(w), the transfer 
function of the retiming circuit filter. 
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Proof: 
Let 
f(j) = elleotte(t)] (69) 


be the input of the bandpass filter and e(t) the input jitter. We 
decompose e(t) into its low-frequency and high-frequency parts, e;(t) 
and ey(t), respectively. We define e,(t) as the component of e(t) with 
frequencies much less than the bandwidth of the bandpass filter. Thus 


—— e(t) = er(t) + en(t). (70) 
Then if | ex(t) | <1 we can write f(t) as 
f(t) = ele] + jen(t) — +--+], (71) 


which is a narrowband process centered at wo. Let Hz(w) be the 
normalized low-pass equivalent of H(w), where 


H,(w) ra H(w + wo) e7 Joo) w > —wWp. (72) 
A(wo) 
Then the output of the filter is* 
g(t) = f(t)*h(t) = e*A(wo)ei {hy (t) xe}, (73) 


where « indicates convolution and h(t), h,(t) are the impulse responses 
corresponding to H(w) and H,(w). The output phase is the phase of 
the term h,(t)*exp(je(t)), plus the static phase shift of e/*, 

Using that [see (72) ] 


i h,(t)dt = H,(0) =] 
and because of the definition of e,(t), we can assume e;(t) to be 
constant compared with h,(t), i.e., 
hy (t)*e* ~ H,(0)e% = eel) 
We therefore obtain 
hy (t)*e? = ef (1 + jey(t)*hr(t)} 

~ eller()+en(O Phy (t) 

a evle(t)*hz(t)], (74) 
Thus, since e(t) is real, 

g (t) = A(wo) ellvot+o(oo)] pile(t)*hi(t)) 


=A (we CEMA] pileott dlwo)+el)"Rehp (ty (75) 
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where 
A(wo)e*"™" ig an amplitude modulation term. Thus the phase 
transfer function is the Fourier transform of w(t) = Re{h,(t)}, that is, 


_ Hiv) + Hi(-w) 


W(o) 5 (76) 
and using (72) we obtain, finally, (68). 
To derive (9) we have made two assumptions: 
(a) Jex(t)| <1 and |ley(t)*h,(t)| «1 (77) 
(b) | ez(t)| = constant. (78) 


When W(w) is going to be used to examine the accumulation of jitter 
in achain of repeaters, the validity of (10) and (11) must be questioned 
for every repeater in the chain. For the Nth repeater, e,(t) represents 
the accumulated jitter appearing in the clock of the (N — 1)-th repeater, 
while e,(t) represents the additional jitter generated by the Nth 
repeater section, and it is always very small. Thus, assumptions (a) 
are valid in general, while assumption (b) is true only if jitter peaking 
does not occur. As is known, jitter peaking occurs if max| W(w) | > 1, 
and in such a case the accumulated jitter grows exponentially. Thus 
relation (68) can be used to evaluate jitter accumulation when we have 
no jitter peaking. If jitter peaking is present, then 


er(t) = Apcos wt, 
where w, is the peaking frequency, i.e., 
| W(w,)| = max| W(w) | 


and (68) can be used only if | e(t) | « 1. This limits the applicability 
of the formula to short chains of regenerators. It is our feeling that in 
the case of jitter peaking, the linear model in (68) will overestimate 
jitter accumulation for long chains, because the nonlinear model will 
shift energy from the peaking frequency band to other bands. Prelim- 
inary simulations appear to agree with the above statement. 

Relation (68) also suggests that 

1. Filters with symmetric ripples in their passband are undesirable 
because they will always create jitter peaking. This is due to the 
normalization factor 2A(wo) in (68). 

2. Even with a monotonic filter, jitter peaking can occur if the data 
frequency, fo, is placed away from the filter’s center frequency f,.’ 
Define the detuning parameter as follows: 

Wo —~ We 
Oa 


B 
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where B is the bandwidth of the low-pass filter H,(w). Then jitter 
peaking usually occurs if a > 1. 

These results are illustrated in Fig. 14, where the phase transfer 
function of a second-order Butterworth is plotted for various a. 


B.2 Timing circuit with prefilter and squarer 


In case a prefilter and a squarer are used we can modify the above 
analysis and also take the above circuits into consideration. 

With wo denoting the baud rate of the received data, the component 
that is going to generate the clock is located at wy/2. Therefore, if e(t) 
is the jitter present in the input data, the component of the input 
located at w,/2 can be represented as 


f(t) a elleotte(t)]/2_ (79) 


Notice that e(t) also contains the jitter generated by the repeater. 

Let us also define as p(t) the output of the prefilter when the symbol 
1 is transmitted, and let P(w) = B(w)e/*™ be the Fourier Transform 
of p(t), i.e., 


p(t) fs Pw) = Blw)ei™, (80) 


Then, using the same assumptions we used in deriving relations (74) 
and (75) and defining 


\w(w)|? 





0 0.5 1.0 1.5 2.0 
FREQUENCY NORMALIZED TO FILTER BANDWIDTH (@/8) 


Fig. 14—Jitter transfer function of a second-order Butterworth filter for various 
values of the detuning parameter a. 
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Plw + wo/2) e7F¥(eo/2) 
B (wo/ 2) 


as the low-pass equivalent of P(w), we obtain 
g(t) = [f(e)*p(t)P*h(e) 
= B(wo /2)eToott2¥ (oo/2)] ellett)"PL()]* p(t) 
= A(c)B(wo/2)eiott2¥o0/2)+ 6(e0)] gile(o*PL(t)*h1 (0) (82) 
Reasoning as in (76) we obtain 


_ Px(w)Hi(w) + Pi(-#)Hi(-e) 


Pi (w) = (81) 


Ww) ao (83) 
Since P(w) is much wider than H;(w), we can assume that 
P,(w)H1(w) = Hz(w), (84) 


which implies that the presence of a prefilter will not change signifi- 
cantly the jitter transfer function of the repeater. However, the jitter 
generated within the repeater may significantly depend on the prefilter 
and squarer. 


AUTHOR 


Christodoulos Chamzas, Diploma degree (Electrical and Mechanical Engi- 
neering), the National Technical University of Greece, Athens, Greece, 1974; 
M.S., 1975, Ph.D., 1979 (Electrical Engineering), Polytechnic Institute of New 
York, Farmingdale, N.Y.; AT&T Bell Laboratories, 1982—. From 1979 to 
1982 Mr. Chamzas was an Assistant Professor with the Department of Elec- 
trical Engineering at Polytechnic Institute of New York, where he is currently 
a part-time Visiting Professor for the Imaging Institute. At AT&T Bell 
Laboratories he has been working on problems in high-speed pseudorandom 
noise generators, the ‘Submarine Lightwave System and adaptive echo cance- 
lers. His primary interests are in signal processing and communication sys- 
tems. Member, Technical Chamber of Greece, Sigma Xi. 


76 TECHNICAL JOURNAL, JANUARY 1985 


AT&T Technical Journal 
Vol. 64, No. 1, January 1985 
Printed in U.S.A. 


Nonlocal Input-Output Expansions 


By |. W. SANDBERG* 
(Manuscript received July 16, 1984) 


In a recent series of papers by this writer on the existence, determination, 
and properties of power-series-like expansions for expressing a nonlinear 
system’s outputs in terms of its inputs, the emphasis is primarily on locally 
convergent expansions. Here we report on related general results concerning 
nonlocal expansions, including in particular material concerning the size of 
the region of convergence. One of the results given provides useful necessary 
and sufficient conditions under which f~' has a generalized power-series 
expansion, where f is a certain important general type of invertible map (that, 
for example, might take one set of complex-valued signals defined on [0, ©) 
into another). 


I. INTRODUCTION 


In a recent series of papers including Refs. 1 through 3 on the 
existence, determination, and properties of power-series-like expan- 
sions for expressing a nonlinear system’s outputs in terms of its inputs, 
the emphasis is primarily on locally convergent expansions. Here, in 
Section II, we report on related general results concerning nonlocal 
expansions, including in particular material concerning the size of the 
region of convergence. 

More specifically, Theorem 1 in Section II gives necessary and 
sufficient conditions under which f~’ has a generalized power series 
expansion (see Section II for the details) when f is an invertible locally 
Lipschitz map between certain general subsets of two complex Banach 


* AT&T Bell Laboratories. 


Copyright © 1985 AT&T. Photo reproduction for noncommercial use is permitted with- 
out payment of royalty provided that each reproduction is done without alteration and 
that the Journal reference and copyright notice are included on the first page. The title 
and abstract, but no other portions, of this paper may be copied or distributed royalty 
free by computer-based and other information-service systems without further permis- 
sion. Permission to reproduce or republish any other portion of this paper must be 
obtained from the Editor. 


77 


spaces. Theorem 2 provides an algorithm for obtaining the expansion 
whenever it exists. 

In Section 2.4 we use Theorems 1 and 2 to prove results concerning 
a system model considered in Ref. 3 and in earlier papers (e.g., Ref. 
2). This model is characterized by five operators: a nonlinear operator 
N and linear operators, A, B, C, and D. The system input v and 
corresponding output w are related by the equations 


y = Nx (1) 
x = Av + Cy (2) 
w = Du + By, (3) 


where x and y, respectively, can be interpreted as the input and output 
of the nonlinear portion of the system.* In Ref. 3 it is assumed that y, 
x, Uv, and w are n-vector valued and defined for t = 0, where as usual n 
is an arbitrary positive integer; here these quantities are allowed to 
belong to a general complex Banach space &, but we shall be interested 
mainly in the case where @ is the space L.. of bounded functions 
considered in Ref. 3. Theorem 4 in Section 2.4 shows that, in an 
interesting and important setting, w can be expressed as a certain 
power series in v that converges for Av € V, where V is any open ball 
in @ centered at the origin such that there is an open subset Vo of @ 
for which (J — CN) (J is the identity operator on & ) is a homeomorph- 
ism of Vy onto V, with (J — CN)" locally Lipschitz on V in the sense 
of Section 2.1. Theorem 4 is actually somewhat more general than is 
indicated above, and it together with its relation to earlier work is 
discussed in Section 2.4.2. A pertinent example is given in Appendix 
B. 


Il. HOMOGENEOUS POLYNOMIALS AND EXPANSIONS FOR MAPS 
BETWEEN SUBSETS OF COMPLEX BANACH SPACES 


2.1 Preliminaries 


Throughout the paper, @ and @ are Banach spaces such that 
there is a linear homeomorphism g of 4 onto &p. (Thus, @ and @ 
are taken to be isomorphic. We shall be interested mainly in the case 
in which in fact @) = & ) The same symbols || - || and 6, respectively, 
are used to denote the norm and zero element in &p as well as in & 
Unless stated otherwise, @) and @ are assumed to be over the field 
of complex scalars. 


* The equations of a very large class of systems can be written in the form (1) through 
(3), with N memoryless. It is not necessary that the equations of the system to be 
studied be given at the outset in the form (1) through (3) (see Ref. 3, Appendices I and 
II). In Ref. 3, y, x, and w are taken to belong to an extended space Lae. For the purposes 
of this paper extended space concepts are not needed. 
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A map g from a nonempty open subset U of @ or of Go into &p or 
@, respectively, is locally Lipschitz on U if for each a € U there are a 
positive number c, and an open ball 6, C U centered at a such that 
lg(u1) — g(ue) |] S cl] uy — uel] for u, and uz in B,. We use d”g(a) to 
denote the mth-order Fréchet derivative (Ref. 4, pp. 149, 181) of g at 
a point a € U, assuming it exists. A sufficient condition for g to be 
locally Lipschitz on U is that it be continuously Fréchet differentiable 
in U. 

For Vo and V any nonempty open subsets of &) and &, respectively, 
H(Vo, V) denotes the set of all homeomorphisms h of Vo onto V such 
that h is locally Lipschitz on Vo. Our results are concerned with this 
class of maps with V assumed to be a c-star about some p € & by 
which is meant that V= {x € @:x =p + W}, where W is a subset of 
@ with the property that zw € W for w € W and any complex scalar 
z with |z| < 1. The concept of a c-star is of importance in studies of 
the region of convergence of abstract power series expansions (see 
Ref. 5, Theorems 26.5.9 and 26.6.1). In this paper we could have 
restricted attention to the case where V is an open ball centered at p, 
which clearly is a c-star about p, but this was not done because no 
significant simplification results. 

Given any positive integer m, by an m-linear map Q from @” into 
#o, we mean that Q(h;, ---, hm) is linear (i.e., additive and homoge- 
neous) separately in each h;. Such a map is symmetric if Q(hi, «++, Am) 
is symmetric in the variables h, ---, hm. A map M(-) from @ into 
&o is called a homogeneous polynomial of degree m if there exists an 
m-linear map Q from @” into @> such that M(h) = Q(h, ---, h) for 
all h*. We now come to a definition of central importance in our 


results. 
For p€ & and VC & an open c-star about p, let “(p, V) denote 


the set of all maps g from V into &p such that there are homogeneous 
polynomials g,,(p, -) of degree m (m = 1, 2, ---), from @ into @o, 
with the properties that 


y &m(p, w) 
m=1 
converges in & for each (p + w) € V, and 


g(p + w) = g(p) + yi E&m(p,w), (pt+tw)e V. (4) 


*This definition of a homogeneous polynomial M(-) is not the same as, 
but is equivalent to, the one given in Ref. 5, that M(zh) = z"“M(h) and M(p + zh) = 
2, Mj(p, h)z/ for p and h in # and any complex scalar z, where the M; do not depend 
on z. Also, the first definition given above describes the same class of maps M(-) if 
“m-linear” is replaced with “symmetric m-linear”. 
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The set #(p, V) is of course a set of maps g that admit a generalized 
power series expansion in the sense indicated. The expansion (4) for 
any g © #(p, V) is unique in the sense that if 


g(p + w) = g(p) + x hn{p,w), (p+tw)EeVv (5) 


[by which is meant, in particular that the sum in (5) converges] with 
each h,,( p, -) a homogeneous polynomial of degree m, then gn(p, -) = 
hn(p, -) for all m. This follows from a simple argument due to Graves 
(Ref. 6, p. 174). (See also Ref. 1, Section 2.7.) 

Finally, we say that g belongs to #;(p, V) if g belongs to #(p, V) 
and for each m there is a continuous symmetric m-linear Q,, from @” 
into &p that depends on p such that g,,(p, h) = Qn(h, --+, h) for all 
h. In particular, then each g,,(p, -) is bounded in the sense that there 
is a positive constant p,, such that ||gn(p, h)|| < pm|l|hl|” for all m and 
h, with p fixed, and every gm(p, -) is Fréchet differentiable on & 


2.2 Inverses of maps in H(Vo, V) and generalized power series expansions 


Our first result, Theorem 1 below, provides a complete characteri- 

zation of those f’s in H(Vo, V), with V a c-star, for which f-' has a 
generalized power series expansion of the type described in the pre- 
ceding section. 
Theorem 1: Let Vo and V be nonempty open subsets of @ and &, 
respectively, with V a c-star about some point p in &. Let f € H( Vo, V). 
Then, f-' € P(p, V) if and only if f is Fréchet differentiable on Vo and 
f~’ is locally Lipschitz on V. In addition, we have f"! © Pr(p, V) 
whenever f-' € P(p, V). 


2.2.1 Proof of Theorem 1 


We first prove the following lemma. 

Lemma 1: Let Vo and V be nonempty open subsets of B and &, 
respectively, with By and & over the same field, either real or complex. 
Let h be an invertible Fréchet continuously differentiable map of Vo 
onto V, with h™ locally Lipschitz on V. Then d(h)(-) (the Fréchet 
derivative of h~') exists and is continuous throughout V. 

Proof of Lemma 1: With q the homeomorphism mentioned in Section 
2.1, define s on Vo by s(x) = gh(x) for x € Vp. It is not difficult to 
verify that s is a continuously Fréchet differentiable invertible map of 
Vo C @p onto the open subset* q(V) of &, and that s~’ is locally 
Lipschitz on q(V). Since the inverses h~' and s~ of h and s, respec- 


* The set q(V) is open by the open mapping theorem, or because it is simply the 
inverse image of the open set V under the continuous map q7}. 
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tively, are related by h7! = s~'q, it follows that the lemma holds in 
general if it holds for @ = Gp. We therefore assume throughout the 
remainder of the proof that @ and &p are the same spaces. 

Let a € Vp be arbitrary, and let chia) and Bra), respectively, be a 
positive constant and an open ball in V centered at h(a) such that 
|h7(u1) — h7"(uz)|| < cael: — uel] for ur and uz in Biya). Let S = 
{x € @:(x + a) € Vo}. Let B denote the open ball in & centered at 
the origin, with the same radius as Baya). 

Define H,:S — @ by H,(x) = h(x + a) — h(a), x € S. The set S is 
open in & Thus, by the continuity of H,, the inverse image H71({) is 
open in S and hence in @&. [The fact that H7z'(8) is an open subset of 
@ is used in connection with (6) and (7) below, and at the end of the 
proof. 

Let h,:Hz'(8) — B be given by h,(x) = H,(x), x © Hz'(8). By the 
invertibility of h, ha is an invertible map of Hz'(8) onto 6. Let hz? 
denote its inverse, and define g:8 — & by 


g(u) = u — dh(a)hz1u + wv, uEB 


for w € 8, where dh(a) is the Fréchet derivative of h at a. 
By the mean value theorem in Ref. 4, p. 160, and the continuity of 
dh(-), 


| ho(2x) — haze) ~ dh(a)(er — 22) 


0 6 
ial 6) 


as max(||2;||, |]z2||) —> 0 with z, ¥ z2.* Therefore, with o > 0 such that 
OCh(a) < 1, there is a 6 > 0 for which {z € @: ||z|] < 6} C Hz1(8) and 


I|halzi) — Ral22) — dh(a)(z: — 22) Il F 
|Z1 = 22|| 


(7) 


for z,; # Ze, ||z:|| < o, and |}z|| < 6. Since ||z]] < caya|lul| for u € B 
and z = h7z'u [because then h(a + z) = u + h(a), which gives z = 
h7 (u + h(a)) — h7*(h(a))], we have 


Ju: — U2 — dh(a)[ho'u, — hous] = 
|hotuy, — ho*uell a 


for u; # Us, || || < p, and ||u2|| < p, where p = min(p¢/a, 5(cna))*), pe 
is the radius of 8, and a is any number in (1, ©). 

Now let u; and uz belong to 6, and define x, and x2 by x; = hztw 
and x. = hz'up. Clearly, h(x; + a) = h(a) + u, and h(x_. + a) = h(a) + 
U2, which gives ||x,; — x2|| < cy a)|| U1 — Ue||. This, together with (8) and 


(8) 


* With regard to Ref. 4, p. 160, ha(z1) — ha(z2) — dh(a)(z; — 22) = [h(a + 21) — dh(a)z] 
— [h(a + 22) — dh(a)zz]. 
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6Ch(a) < 1, shows that the map g defined above is a contraction on the 
set So = {u€ FB: |\ul| < p}. Since hz'(6) = 8, it is easy to see that for 
all w € F with ||w|| S po and po > 0 sufficiently small, g maps Sp into 
itself. By the contraction-mapping fixed point theorem, g has a unique 
fixed point* in So, and thus there is a unique solution u € So of 
dh(a)hz'u = w, for each such w. 

By the linearity of dh(a), we see that for each w € @ there is a 
unique x © @ such that dh(a)x = w, the uniqueness following from 
the fact that h, maps an open ball in @ centered at # into So [recall 
that h,(6) = 6, that h, is continuous, and that H71(8) is open in @].* 
This shows that dh(a) is an invertible map of @ onto &. By the 
boundedness of dh(a), dh(a)~ is bounded (Ref. 9, p. 119). Since a € 
V. is arbitrary, and dh(-) is assumed to be continuous, it follows (Ref. 
4, p. 273) that h~! is continuously Fréchet differentiable on V, which 
proves the lemma. 

Returning now to the proof of the theorem, suppose initially that 
f-! € P(p, V). Since the series associated with f~’ is a (G)-power 
series in the sense of Ref. 5, p. 773, this series is (G)-differentiable 
(i.e., Gateaux differentiable) in V (Ref. 5, p. 773). By the (G)-differ- 
entiability and continuity of f-' (which implies that f~! is analytic), it 
follows Ref. 5, Theorem 3.17.1 that f~! is Fréchet differentiable 
throughout V. Using the following lemma,’ which is proved in Appen- 
dix A and which is used also in Section 2.3, d(f~')(-) is continuous. 

Lemma 2: If g is a Fréchet differentiable map of a nonempty open 
subset V of & into a complex Banach space &,, then g is twice Fréchet 
differentiable on V. 


* This part of the proof is along the lines of proofs of the classical implicit function 
theorem (see, for example, Ref. 7, pp. 194-5). Lemma 1 is related to a theorem stated 
in Ref. 8, p. 165, but the argument given there does not prove the theorem. (In the 
terminology of Ref. 8, it is not shown that the definition of a linearization leads to the 
limit given. However, the proof above shows that the theorem in Ref. 8, p. 165 is true 
under the additional hypotheses that L, exists and is continuous in a neighborhood of 
xX = X.) 

* A modification of an argument given in Ref. 8, p. 165 could also have been used to 
obtain the uniqueness. 

* With regard to Lemma 1 and the condition in Section 2.1 that there is a linear 
homeomorphism g that maps 4 onto &p, we note that, in the absence of that condition, 
the existence of dh(-) and dh~'(-) as indicated in the lemma implies (see Proposition 1 
in Section 2.3.1) that the condition holds. This shows that there is no loss of generality 
in assuming that the condition is met. 

§ This lemma, which is an analog of a standard proposition in the classical theory of 
functions of a complex variable, is probably a known result, but we have not encountered 
it in the literature. However, for related material concerning Gdteaux variations, see 
Ref. 5. The lemma provides additional understanding concerning some results proved 
in Ref. 1 and other recent papers. While it does not strengthen these results, it shows 
that some hypotheses actually follow from others that were introduced (see, e.g., Ref. 1, 
Theorem 2, which is an early result concerning nonlocal expansions). 
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Thus f~' is locally Lipschitz on V, and, by Lemma 1, df(-) exists on 
Vo. Using the analyticity of f~’, by Ref. 5, Theorems 3.16.2, 26.3.4 and 
26.3.6, one has f-' € Pp(p, V). 

Now assume that f is (F)-differentiable, and thus continuously 
(F)-differentiable, on Vo and that f~’ is locally Lipschitz on V. 
Using Lemma 1, d(f~')(-) exists in V. In particular, f~' is (G)-differ- 
entiable in V, and by Ref. 5, Theorems 3.16.2 and 26.3.5, we have f-' € 
@\(p, V), which completes the proof of the theorem. 


2.2.2 Comments 


Proposition 1 in Section 2.3.1 and the proof of Theorem 1 show that 
under the hypotheses of the theorem, f-' € #(p, V) implies also that 
df(a) is a homeomorphism of @» onto @ for every a € Vo.* 

For the extreme case in which both & and & are just the space of 
complex numbers with the absolute value norm, the condition that V 
is an open c-star about p reduces to the requirement that V is an open 
disk in the complex plane centered at p. In that case, by standard 
results in the classical theory of functions of a complex variable, f~' € 
Pp, V) implies that the ordinary derivative (f~')’(z) exists for 
z € V. In addition, since f-' is one-to-one on V, it follows from a 
known result (Ref. 10, Theorem 16-23) that (f~')’(z) # 0 for z € V, 
which shows that f’(-) exists on Vo. Similarly, using Ref. 10, Theorem 
16-23, it follows from the hypothesis that f’(z) exists for z € Vo that 
f-1 € A(p, V).* Theorem 1 can be viewed as a Banach space relative 
of these propositions. 


2.3 Construction of the expansion of f~' 
Here we give an algorithm, along the lines of Theorem 2 of Ref. 1, 
for expanding f'. 


Theorem 2: Let the hypotheses of Theorem 1 be met, and assume that 
f-! © Alp, V). Then for each | = 1, 2, ---, the lth-order Fréchet 
derivative d'f(x) exists for x € Vo, and 


fp +h) =f"(p) + X Sm(p,h), (p+h)EV (9) 


where 


&i(p, h) = (df[f"(p)]) th 


* And this adds to material in Ref. 1, Section 2.6 concerning the necessity of a certain 
invertibility condition. 

tIn these observations, the hypothesis that f is locally Lipschitz is not needed. On 
the other hand, the Lipschitz condition is obviously a consequence of f-' € #(p, V). 
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(the inverse exists by Proposition 1, below), and 
&m(p, h) = —(df[f"(p)))* p> (1!) 


-> dif (p) gs,(p, h)gp,(p, h)-+-82,(p, 2), 


Ry-thyt-+-+ky=m 
j 
m22.* (10) 
2.3.1 Proof of Theorem 2 


Proof: Since by Theorem 1 and Proposition 1 (which appears below), 
f- € P(p, V) implies that df(a) and [df(a)]~’ exist for each a € Vo, 
Lemma 2, the proof of Theorem 2 of Ref. 1, and Theorem 3.16.2 of 
Ref. 5, show that Theorem 2 holds. 


Proposition 1: If Vo and V are nonempty open subsets of @) and &, 
respectively, and f is a homeomorphism of Vo onto V such that f and 
f~ are Fréchet differentiable on their respective domains, then df(a) 
is a homeomorphism of &o onto & for any aE Vo. 

The proposition is a well-known result (see, for example, Ref. 
11, p. 175, Problem 6) provable using the relations f~'[f(x)] = x and 
fif-'(y)] = y for x and y in Vp and V, respectively, and the chain rule 
(see Ref. 11, pp. 171-2) for differentiating a composite function. 


2.4 Theorems concerning the system model 


In this section attention is focused on the system model described 
in Section I. We use J to denote the identity map on &. 

Let A, B, C, and D be linear maps of @ into itself, with B and C 
bounded. Let N be a map from a subset S of @ into & for which 
there are nonempty open subsets Vo and V of @& such that Vp C S, V 
is a c-star about some point p € &, and (J — CN) is a homeomorphism 
of Vo onto V, with (I — CN)" locally Lipschitz on V. (It is not difficult 
to give important examples in which these hypotheses on N are met. 
This is illustrated in Appendix B.) 

Let _7 (for “set of inputs”) denote the collection of all v € @ such 
that Av € V, and assume that there is a up € & for which Avo = p. 
We see that for each v € _Z there are unique x, y, and w in Vo, & and 
@, respectively, such that 


x= Av+ Cy 
w = Du + By 
y = Nx. 
*In (10), Ditesonak: -m denotes a sum over all positive integers k;, ---, kl that add 
to m. a 
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Let W be the map from _¥ into @ defined by the condition that 
W(v) = w when v € _Z In other words, let W be given by 


W(v) = Dv + BU — CN)"'Ap, ve Z% (11) 


Since N differentiable on Vo implies (by Lemma 2) that it is 
continuously differentiable on Vo and thus locally Lipschitz there, by 
Theorem 1 and the boundedness of C we have the following: 

Theorem 3: Let N be Fréchet differentiable on Vo. Then (I - CN)": € 
Gp (DP, V). 

Theorem 3 shows that, under merely the condition indicated, W has 

an expansion about vp in terms of homogeneous polynomials that is 
valid* for all v € _Z With regard to actually determining the expansion 
for W, we have the following. 
Theorem 4: Let N be Fréchet differentiable on Vo. Then for each | = 2, 
3, --- the Fréchet derivative d' N(-) exists in some open neighborhood 
% of the point xo in Vo that satisfies x» — CN(xo) = Avo, and for 
v € _Z we have 


Wi(v) = BN(xo) + Du + y BY (m)(V — Uo) (12) 
m=1 


where Yq), Yi), --- are the homogeneous polynomials defined by the 
relations* 
X(y(v — vo) = [I — CdN (xo) 'A(v — vo), 


Xim(v ~ ve) =F - CANGOTC YY" 


y d'N(xo)X(z,(V — Vo) +++ X~@(V— vo) (13) 


Rytkot---+hj=m 
k;>0 
for m = 2, and 


Yim (V — Vo) = ¥ (Il!) »> d'N(xo)Xx,)(V — Vo) 
l=] Ryt+kot+---+k=m 


kj;>0 


oe - Xn) (v Uo) (14) 
form2 1. 


2.4.1 Proof of Theorem 4 


Proof: Let @? denote the Banach space @ X @ with norm max 
(l-], l-I). Define G: Vo x @ > @? by 


* Here the boundedness of B is used to ensure that )%-0 Bgm[p, A(v — vo)] converges 
to BYn-0 8mlp, A(v — Vo), where Yn=0 &m[p, A(v — vo)] is the series for (I ~ CN)" Av. 
t The inverse of J — CdN(xo) exists; see Section 2.2.2. 
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G\(Pi, P2) = pi — CNp; 


Go(pi, De) = Npi — Pe 


for (pi, D2) € Vo X & With 8 a nonempty open ballin @ centered at 
6, let 


g = {( pi, P2) E Vo x B:G(pi, Pe) = (q1; q2), (q1; 2) € V x B} 


fie., let D = G(V x B)]. We see that F is open in GB? (because V 
and @ are open, and G is continuous and defined on an open subset of 
@*), and that V x Bis a c-star in Z? centered at (p, 6). The restriction 
F of Gto @ is a differentiable homeomorphism of J onto V x 8. 
Using this fact, it is not difficult to show that a proof of Theorem 4 
can be obtained by proceeding as in Ref. 2, proof of part (ii) of Theorem 
1, but with our Theorem 2 with @ = &p employed instead of Lemma 
1 in Ref. 2.* (It is easy to verify that for F as described above, F~’ is 
locally Lipschitz.) 


2.4.2 Discussion 


Under the conditions on N of Theorem 4, N has the representation 
N(x) = N(x0) + & (I!)~td'N(xo)(x — xo)’ (15) 
i=l 


in which the series converges (for example) uniformly for x in some 
sufficiently small ball in S centered at xo,‘ and of course (15) provides 
an important interpretation of the maps d'N(xo) that appear in (13) 
and (14). 

It is often useful to observe that the operator (J — CN)~ in (11) can 
naturally be identified with the “feedback part” of the system repre- 
sented by (1), (2), and (3) (see Ref. 3, Figs. 1 and 2). Aside from 
considerations concerning the differentiability hypothesis on N, and: 
the existence of a Up as described, Theorem 4 shows that the series 
representation given by (12) holds for Au € V, whenever V is an open 
c-star about some p such that the equation x — CNx = u of the feedback 
portion is, so to speak, uniquely locally-Lipschitz-solvable in some 
open subset of S for every u in V. Implicit in this is the assumption 
that @ is a complex Banach space. Since the inputs and outputs of 
most nonlinear systems of direct interest are real valued functions, 
the main point of Theorem 4 with regard to applications is that an 


* See also Ref. 3, Theorem 1, with regard to the expression for the Yn. 

* See Ref. 5, Theorem 3.17.1, and Ref. 11, p. 198. Also, notice that the convergence 
to N(x) of the right side of (15) for x» = p and x € V is implicit in Theorem 4, since C 
and D can be taken to be the zero map and A and B can be assumed to be the identity 
map. 
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expansion exists, and can be constructed as specified, when the original 
system equations can be “complexified” so that the hypotheses con- 
cerning Vo, V, and (I — CN), and of the theorem, are met. It is not 
difficult to give specific examples. One example is provided by the 
material in Appendix B. 

Theorem 4 with N(6@) = p = Up = Xo = 6 bears directly on the main 
result in Ref. 3, which concerns expansions involving iterated inte- 
grals, for the case in which @ is the set L.. of bounded complex n- 
vector valued functions on the interval [0, 0).°* Our theorem shows 
that the locally convergent expansions described there converge in 
fact for Av € V for any V and Vo such that our solvability hypotheses 
with p = 6 are met, assuming merely that S = I, where I is the domain 
of definition of N in Hypothesis B.2 of Ref. 3. 

In this connection, if @ = L. and N is memoryless in the sense 
that (Ns)(t) = n[s(t), t], ¢ = 0 for s € S, with S = [ and 7 and T, 
respectively, a function and a domain of the kind described in Ref. 3, 
Hypothesis B.2, and if, for simplicity, n = 1, then [d'N(xo)X«,) 
(v — Uo) +--+ X(ay(V — Uo)](t) is just the product M1(t)[X,)(v — vo)](t) 
-++[X¢)(V — vo)(t)] for t = 0, where M1(t) = d'n(z, t)/d2'| 2-4) (see 
Ref. 2, Lemma 3). 
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APPENDIX A 
Proof of Lemma 2 


Proof: Choose a € V and, using the continuity of g, let r and M be 
positive numbers such that Z C V and ||g(x)|| < M for x € G where 
D= {x € @:|\x — all <r}. Since g is analytic on ZF in the sense of 
Ref. 5, Definition 3.17.2, there are continuous symmetric m-linear 
maps Q,, from 4” into &,, which depend on a, such that 


g(a) + Y (m!)"Qm(k, «++, R) 
m=1 
converges absolutely to g(a + k) for ||k|| <r (see Ref. 5, Theorems 


3.16.2, 26.3.4, 26.3.6, and 26.6.6). 
With h € @ consider the formal series 


ims 


(m!)~"'mQn(h, sg h, -) (16) 
1 


whose terms are bounded linear maps from & to &,. Choose any 
GB, > 1 and let 8 = (8,/r). Using an analog (Ref. 5, Theorem 3.16.3) of 
the Cauchy bounds and the m-linearity of Qn, 


WQm(k, +++, R)l| < M(m!)(B RII)” 
for k © &. Since (see Ref. 16, Proposition 2.2.11) 
sup{ || Qm(Fi, «++, Rm){l = Will <1, +++, [Anil <2) 
<m™(m!)“'sup{]Qn(R, ++, Rl: WRI <1) (17) 


we have 
Qm(Ri, +++, Rm)|| < m"™MB™ [Ril] --+ Rm (18) 
for any ki, ---, Rk, in & By Sterling’s formula for m!, 
m! > (22) /?2m'/?2m™e—™. (19) 


Thus, by (18) and (19) one has for the mth term in (16), 
lm!) "mQn(h, «++, h, Dll 
= sup{||(m!)"'mQ,,(h, «++, h, RD: WRI <j 
< (20)-'?m'?M (eB) |\eBh||"-”. 
It easily follows that there is a positive r; < r such that (16) converges 


absolutely to a bounded linear map L(a, h) for ||h|| <r:. Since L(a, h) 
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is an (F)-power series in h in the sense of Ref. 5, Theorem 26.64, the 
map L(a, -) is analytic and hence Fréchet differentiable (Ref. 5, 
Theorem 3.17.1) on {x € @: ||x|| <1r,}. In particular, L(a, -) is Fréchet 
differentiable at 6. 

In addition, for ||h|| < ar; and ||6|| < ar, where a € (0, 1/2), 


g(a +h + 6) — g(at+ h) — L(a, h)é| 
S d, ll (mm!) [Qn (h a 6, te A+ 6) > Qn (h, wre os h) 


po mQn(h, ery h, 6) Il 


IN 
ims 


: (m!)772” || Qm{|(ar)" [5], (20) 
where || Q,|| is the left side of (17). Using (18) and (19), it is a simple 
matter to verify that for a sufficiently small the extreme right side of 
(20) is o(||6]}) as ||6|] — 0. This shows that dg(a + h) = L(a, h) for all 
h in some neighborhood of 8, and thus that dg(-) is Frechet differen- 
tiable at the arbitrary point a € V. 


APPENDIX B 
An Example 


Here we give a simple example of a map (J — CN) and associated 
sets Vo and V that meet the conditions at the beginning of Section 
2.4. The construction of the example involves a contraction-mapping 
technique that is well known and often useful. We shall use L.. to 
denote the set of bounded complex-valued Lebesgue measurable func- 
tions defined on [0, ©) with the usual sup norm.* 

Let S = @ = L., and let us take C to be an operator, such as a 
convolution operator, with induced norm p. We assume that p > 0. Let 
N be defined by (Nx)(t) = x(t)? for t = 0 and x € @& and take F: 
B&B — & to be given by Fx = u + CNx for every x € & where u is an 
element of &. 

We note that for xo in the closed ball @(r) of radius rin @ centered 
at 6, dN(xo)h(t) = 3xo(t)*A(t), t = 0, for arbitrary h © &, and therefore 
|| CAN (xo) || < 3pr?. Thus, F is a contraction mapping on @(r) for 0 < 
r<ro, where ro = (3p). Also, F takes @(r) into itself if 


[ul] + sup{|CNx]:x € B(r)} <r, 


*In other words, here L. is the set L. in Ref. 3 with n = 1. (The definition of L. 
ordinarily found in the system-theoretic literature involves instead the essential sup 
norm. That our L. in this paper, or in Ref. 3, is complete follows from the Cauchy 
criterion for uniform convergence and the standard proposition that pointwise limits of 
measurable functions are measurable.) 


NONLOCAL EXPANSIONS 89 


which is met when ||u|| < g(r), where g(r) = r— pr®. These observations 
and the monotonicity of g on (0, ro) motivate the following. 

Take V to be the open ball in @ of radius q(ro) centered at 9. Given 
any u € V there is an r; € (0, ro) such that ||u|| < q(r) for all rE 
[71, To). Hence, by the contraction-mapping theorem, for any such u 
and all associated r there is a unique solution x € G(r) of the equation 
x — CNx = u. Also, for any open ball 8, in V centered at some point 
a, there is an r € (0, ro) such that ||u|| < q(r) for u € §,, and it easily 
follows that the solutions x, and x2 corresponding to any u, and uz in 
Bq satisfy |x, — x2|| < (1 — 3pr?)7? |] uw, — uel. 

Thus, with Vo the intersection of the open ball 8o(ro) of radius ro 
centered at @ with the inverse image of V under (I — CN), we see that 
V> is open and that (J — CN) is a homeomorphism of Vo onto the open 
c-star V about the point 6, with (I — CN)” locally Lipschitz on V.* 
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High-speed (2 4.8 kb/s) echo-cancellation-based full-duplex direct distance 
dialing modems usually have to deal with two echos: the near echo, which is 
generated at the modem location, and the far, or talker, echo, which has been 
looped back to the modem after passing through a carrier system. The near 
echo propagates through a channel that is essentially linear, and thus it can, 
at least in theory, be perfectly canceled by an echo canceler. On the other 
hand, the far-echo channel is generally plagued by impairments that can 
seriously degrade the performance of an echo canceler. In this paper we study 
the effects of these channel impairments on the performance of an in-band 
data-driven echo canceler. This echo canceler has been found to be particularly 
well suited for full-duplex voice-grade data transmission applications. Both 
analytical and real-time experimental results are presented. It is shown that 
frequency offset, even in small amounts, is by far the most damaging of the 
channel impairments that are commonly encountered in carrier systems. The 
degradation of performance due to phase jitter can be significant. However, 
this can only happen under simultaneous worst-case conditions of phase jitter 
and signal power levels, and these cases might not be statistically significant. 
Worst-case nonlinearities in the echo channel do not degrade substantially 
the performance of the echo canceler. 


I. INTRODUCTION AND SUMMARY 


In this paper we study the effects of channel impairments on the 
performance of a data-driven in-band echo canceler. This echo can- 
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celer is intended to be used in two-wire, high-speed (> 4.8 kb/s), full- 
duplex, voice-grade data transmission. Details about this application 
for 4800 b/s Direct Distance Dialing (DDD) operation are given in 
Ref. 1. Performance degradations due to finite precision effects in the 
digital implementation of the echo canceler are studied in a companion 
paper.” 

The use of echo cancellation in two-wire, full-duplex data transmis- 
sion will be described with reference to Fig. 1.6 Figure la shows a 
simplified version of a typical connection over the switched network. 
Echoes arise because of impedance mismatches in the hybrid couplers 
that make the connections between two-wire and four-wire transmis- 
sion facilities. Consequently, some energy leaks directly through the 
first hybrid encountered by the transmitted signal. This signal is called 
the near echo. Similarly, some energy leaks through the hybrid at the 
other end of the four-wire circuit and is looped back to the modem 
through the carrier system. This signal is called the talker or far echo. 
These two echoes are added to the signal transmitted by the far-end 
modem, and they appear as interference at the receiver’s input. 

The echo canceler in Fig. 1b synthesizes a replica of the channel 
traversed by the echoes. It can take its input from various points of 
the transmitter. For example, a voice-type canceler takes its input 
from the output of the transmitter, whereas a data-driven canceler 
gets its input directly from the data symbols at the input of the 
transmitter. Under ideal conditions, if the canceler and the echo 
channel have the same inputs, they should also have the same outputs, 











FOUR-WIRE 
CARRIER CIRCUIT 


TWO-WIRE 
LOOP 


(a) 


DATA 
TR NEAR ECHO 
ECHO 
FAR ECHO 


(b) 


Fig. 1—(a) Typical dialed connection. (b) Use of data-driven echo canceler at station 
location. 
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and the signal, after subtraction at point A, should consist only of the 
wanted far-end data signal. These ideal conditions exist, to some 
extent, for the near-end echo because the hybrid introduces mainly 
linear distortion (for which the canceler can compensate perfectly if 
its memory span and the digital precision are large enough). However, 
this is not true for the far echo that has propagated through a channel 
that is generally time varying and nonlinear. The updating algorithms 
and structure used in the basic canceler can only compensate partially, 
if at all, for these impairments, and some residual echo will appear at 
point A in Fig. 1b. Some of the effects of these channel impairments 
for voice-type cancelers are studied in Refs. 7 through 11. 

Because the near echo and the far echo have different characteris- 
tics, we will find it convenient to break the canceler into two parts, a 
near canceler and a far canceler. The requirements for these two 
cancelers are quite different. Due to the delay characteristics of the 
far-echo channel, the far canceler requires a larger memory span than 
the near canceler. However, under worst-case conditions of hybrid 
leakage and signal levels, the near canceler has to provide a much 
larger echo attenuation than the far canceler. In both cases these 
attenuations generally have to be achieved under “double-talking”* 
(full-duplex transmission) conditions and, in the case of the far can- 
celer, in the presence of time-varying and nonlinear echo-channel 
impairments. 

Our purpose in this paper is to study the performance degradation 
of the far canceler in the presence of the impairments that are most 
commonly encountered in carrier systems. The echo canceler used in 
the study is the so-called data-driven, Nyquist, in-band canceler. Its 
structure is derived in the next section. The experimental and analyt- 
ical work presented here consisted of studying the signals at point A, 
in Fig. 1b, when different types of impairments were inserted in the 
far-echo channel. In the experimental setup, the echo channel was 
simulated by using the appropriate laboratory equipment, and real- 
time performance measurements were achieved by using an echo 
canceler implemented on an in-house developed digital signal proces- 
sor. We now summarize our findings. Frequency offset, even in small 
amounts (a fraction of a hertz), is found to be the most damaging of 
all the impairments. Possible corrective actions are studied in a 
forthcoming paper. The performance of the canceler can also be 
significantly degraded by phase jitter. However, in practice, this can 


* The signal at point A in Fig. 1b is used to adapt the echo canceler’s tap coefficients. 
In steady-state operation this signal will consist mostly of the far signal, which will thus 
act as a strong noise component in the adaptation algorithm.’ An arrangement proposed 
recently allows, in principle, the elimination of the double-talker from the adaptation 
algorithm.” 
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only happen under certain simultaneous worst-case (round-trip) con- 
ditions of jitter and signal power levels, and these cases might not be 
statistically significant. Worst-case nonlinearities did not degrade 
substantially the performance of the canceler. 

The paper is organized as follows. The structure of the echo canceler 
used in this study is derived in the next section. Its convergence 
properties in the presence of a double-talker are analyzed in Section 
III. The effects of phase jitter and frequency offset are studied in 
Sections IV and V, respectively. Finally, in Section VI we present real- 
time performance results obtained with an echo canceler implemented 
on a bit-slice processor. 


Il. IN-BAND ECHO CANCELER STRUCTURE 


The echo canceler structure derived in this section is of the Nyquist 
or interpolating type. That is, it cancels the echo at all frequencies. It 
is also called an in-band canceler because it synthesizes passband 
filters rather than equivalent baseband filters. This canceler is partic- 
ularly attractive because it is less complex to implement than either 
the voice-type canceler or the Nyquist data-driven canceler described 
in Ref. 6. 

A two-dimensional (in-phase and quadrature) modulated signal is 
generally represented by the expression 


s(t) = Re 12 Ang(t — nner (1) 


where A, = a, + jb, is the discrete-valued multilevel complex symbol 
to be transmitted, g(t) is a Nyquist pulse, 1/T is the symbol rate, and 
w,/27 is the carrier frequency. In the usual case where the highest 
frequency component in g(t) is smaller than the carrier frequency, the 
complex signal in brackets in (1) is an analytic signal Z(t), where 


Z(t) = s(t) + js(t) = ¥ Ang(t — nT)e, (2) 


and where s(t) is the Hilbert transform of s(t). Equation (2) can be 
rewritten as 


Z(t) = Y Ane’"* g(t — nT el" (3) 
Z(t) = } A,R(t — nT), (4) 
where 
AL = A,e@"T (5) 
R(t) = g(te*“. (6) 


94 TECHNICAL JOURNAL, JANUARY 1985 


When the signal s(t) is transmitted through a channel with impulse 
response h(t), the analytic signal corresponding to the output signal 
is 


Z(t) = Z(t)*h(t) = Y AARi(t — nT), (7) 


where * denotes convolution and 
R(t) = R(t)*h(t), (8) 


and where Z;(t) and R,(t) are analytic signals. The signal at the output 
of the channel is the real part of Z;(t), ie., 


si(t) =  [aan(t — nT) — b,A(t — nT)], (9) 


where R,(t) = r,(t) + jF,(t). 

Thus, the echo s(t) is obtained by feeding the symbols aj, and b;, to 
in-phase and quadrature bandpass filters with impulse responses r;(t) 
and 7,(t), respectively. This suggests the following structure for a 
digitally implemented data-driven echo canceler. We can feed the 
rotated symbols aj, and b/, to two transversal filters with variable tap 
coefficients and tap delay spacings of JT’. The signal s,(t) in (9) is 
sampled at a rate 1/T’, which will be assumed to be at least twice the 
highest frequency of s;(t), so that the wanted far signal, which is added 
to s,(t) in full-duplex operation, can be reconstructed after cancella- 
tion. A standard Mean-Squared Error (MSE) criterion can then be 
used to adapt the tap coefficients. Ideally, after adaptation, the tap 
coefficients of the two transversal filters converge to the sampled 
values of r,(t) and F,(¢) in (9). The echo canceler structure is shown 
in Fig. 2. Notice from (4), (5), and (6) that the same structure can be 
used to implement the transmitter, in which case the tap coefficients 
are fixed.’* 

Referring again to Fig. 2, the transmitter generates the symbols a, 
and b; at the symbol rate 1/7, which is always smaller than the 
sampling rate 1/T’. Therefore the canceler’s delay lines in Fig. 2 will 
only be sparsely filled. In fact, they will contain L — 1 zeros for each 
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Fig. 2—In-band data-driven echo canceler structure. 


CHANNEL IMPAIRMENTS OF ECHO CANCELER 95 


SUBCANCELER (c_,d_g) 





Fig. 3—(a) Modified canceler structure. (b) Subcanceler structure. 


nonzero symbol, where L = T/T’ is the number of sampling periods 
per symbol interval.* A direct implementation of this canceler struc- 
ture would waste computation power since multiplication by zero need 
not be performed. A more efficient, but equivalent, structure can be 
derived by observing that a different subset of coefficients is used for 
the computation of the output when the symbols move down the delay 
line in Fig. 2. There are L such subsets that are each used once per 
symbol interval. All the subsets of coefficients are correlated with the 
same input vectors (a7, an-1, ---) and (b;, bj-1, ---) in a given symbol 
period. Therefore the canceler can be considered as a parallel combi- 
nation of L “subcancelers” having all the same delay line but different 
sets of coefficients, as shown in Fig. 3a. One of the subcancelers is 
shown in Fig. 3b. The taps of the transversal filters are now spaced at 
T rather than T’, and none of the entries in the delay lines are zero. 
The ith subcanceler consists of the tap vectors c; and d; whose elements 


* Theoretically many other interpolation schemes can be envisioned. However, the 
scheme described here leads to the simplest possible practical implementation of the 
echo canceler. 
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are C_n+itke, And G_nsise,, respectively, where 1 = 0 --- L — 1. The 
outputs of the canceler are generated at the sampling rate 1/T’ by 
computing, in a cyclic fashion, the outputs of the L subcancelers shown 
in Fig. 3a. This model applies also to the echo channel and we can 
consider that it consists of a parallel combination of LZ subchannels. 
Convergence of the whole canceler can then be achieved by having 
each subcanceler converge to the corresponding subchannel. The al- 
gorithms used during adaptation are described in the next section. 


II. CONVERGENCE PROPERTIES IN THE PRESENCE OF A DOUBLE- 
TALKER 


In order to adapt the tap coefficients of the canceler, we will 
separately minimize the MSE between the outputs of each subcanceler 
and the corresponding echo subchannel. That is, each subcanceler’s 
output is computed once per symbol interval, an error is derived, and 
a stochastic-gradient algorithm is used to update the tap coefficients. 
Nyquist cancellation is obtained by cyclicly repeating these operations 
at the sampling rate 1/T’ for all the subcancelers. We analyze this 
adaptation scheme by assuming that the subcancelers adapt independ- 
ently. This assumption was found to be in excellent agreement with 
our experimental results. For simplicity of notation, we will also 
assume that the inputs of the echo canceler are the original symbols 
a, and b, rather than the rotated symbols a}, and b7. We now define 
the following quantities, at the nth symbol instant: 


al = [an, An-1, An-2, ---] = in-phase data vector. 

b? = [b,, bp-1, bp-2, ---] = quadrature data vector. 
‘as 

C; = [C_n+i, C-N+itL, C-N+i+2L) °° -] 


= vector of in-phase tap coefficients of ith subcanceler. 
d? = [d_n+i, d_-nsitt, d-nsitar, °° | 


= vector of quadrature tap coefficients of ith subcanceler. 


ri = (n(0:T’), (0 + iT’), r(2T + iT’), ---] 
= vector of in-phase samples of ith echo subchannel. 
ro; = [ro(iT’), ro(T + iT’), ro(2T + iT’), «+ +] 
= vector of quadrature samples of ith echo subchannel. 
The superscript 7T' designates transposed vectors where all these 


vectors are assumed to be infinite. For finite-length cancelers shorter 
than the echo-channel impulse response, we simply insert zeros on the 
right in the definitions of e7 and d}. The outputs, s.(nT + iT’) and 
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s(nT + iT’), of the ith subchannel and the ith subcanceler at time 
nT +1T’ are 


se(nT + iT’) = alr; + bore + Eni (10) 
s(nT + iT’) = aze; + bid, (11) 


where é,,; is an additive interference signal that is uncorrelated with 
the signal to be canceled. This interference will generally consist of 
the desired (far-end) data signal and some additive noise. The error 
e,,; between the outputs of the ith subchannel and the ith subcanceler 
is 


Cn; = S(nT + iT’) — s.(nT + iT’), 1=0,1,---,L-—1, (12) 
and we want to minimize the MSE 
E; = ((eni)”) = ([se(nT + iT’) — s.(nT + iT’)]’) 
= (fan(rii — ¢;) + bara — di) + &,,]*), (13) 


where (-) denotes the expectation of the quantity inside the brackets. 
Note that E; cannot be smaller than the irreducible noise (£3,;). 

The Minimum MSE (MMSE) is achieved when the N; complex tap 
coefficients of the ith subcanceler are equal to the corresponding 
subchannel complex sampled values. The MMSE is given by 


min E;=A ¥ [ri(kT + iT’) + r3(kT + iT’)] + (&?), (14) 
k>N; 


where the data-symbol power is given by A = (a”) = (b’”). 

This MMSBE is obviously not the same for all the subcancelers. The 
adjustment algorithms for the updating of the subcanceler tap coeffi- 
cients are obtained by taking the gradient of the MSE in (13) with 
respect to the tap vectors ¢; and d; 


0E; 
2E = 2 ane) = Bu © (15) 
OE; 
ad, = —2(bneni) = ro — d;. (16) 


As is usual in practice, the gradients, with respect to the squared error 
rather than the MSE, are used for the adjustment of the tap coeffi- 
cients. The corresponding stochastic tap adjustment algorithms are 
then 


Ch+1,i = Cri + AAanen,i (17) 
Gn+1i = d,.i SB abreni, (18) 


where a is the step size of the adjustments andi = 0, 1, ---, L—1. 


98 TECHNICAL JOURNAL, JANUARY 1985 


Inspection of (11), (17), and (18) shows that each subcanceler requires 
a total of about 4N multiplications and additions for the filtering and 
updating operations. Therefore the implementation of the whole can- 
celer requires 4LN multiplications and additions per symbol period. 
An analysis similar to those given in Refs. 4, 6, and 14 can be used to 
study the convergence properties of the MSE as a function of time. 
As is usually the case, the analysis assumes that the data vectors a, 
and b, in the subcancelers are uncorrelated between successive tap 
adjustments. Under these conditions it is shown in Appendix A that 
the MSE decreases as 


(e2;) = (1 — 2aA + 2a?N;A”)"(e6,:) 


1 — (1 — 2aA + 2a°N;A?)” 


2. 
1 — (1 — 2ad + 20°N,A®) AE 8) 


where N; is the number of complex taps in the ith subcanceler. 

In the derivation of (19) it was assumed that the first term on the 
right in (14) was zero. That is, it was assumed that the canceler was 
long enough to cover the memory span of the echo channel, and that 
no degradation was introduced by the finite precision in the digital 
implementation. Effects of finite precision are studied in detail in 
Ref. 2. 

For the expression in (19) to converge, we require 


11 — 2aA + 2a°N;A?| <1 (20) 
so that the step size a has to satisfy 





O<a< (21) 


N;A° 
The step size that provides the fastest speed of convergence to the 
corresponding steady-state MSE is obtained by setting the derivative 
of the expression in (20) to zero, i.e., 


1 


Qopt cad 2N;A ) (22) 


or one-half the maximum step size. From (19) it is clear that the ith 
subcanceler’s steady-state MSE for a given step size is given by 


Pree (i) 
and the overall steady-state MSE averaged over a symbol period is 
ae) 
Bee ee ee : 
(en) =F BT ond’ eo 
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If we assume that all the subcancelers have the same number N of 
complex taps and if we define 


aL ye 25 
(ee oe) (25) 


as the average power of the interfering signal in a symbol period, we 
can rewrite (24) in the form 
2 
Licgethey Se ND 
E = (ea) i oN (26) 
The signal-to-noise ratio (s/n) achievable in front of the receiver 
can be derived from (26). The uncorrelated term (£*) consists of the 
noise in the channel with power P, and the far signal with power Ps. 
The residual MSE E is the sum of Ps and some interfering “noise”, J, 
defined by 


7 _ aNAPs + P, 
obo Pie a (27) 


The steady-state received signal-to-noise ratio is then 


Ps 1- aNA 
eo GNA & PUPS co 
Notice from (22) and (28) that the maximum achievable s/n is 0 dB 
when the optimum step size for speed of convergence is used, and the 
channel noise is assumed to be zero (P,, = 0). Although the noise in 
(22) through (28) is not Gaussian, this expression has proved to be 
very useful for predicting the echo canceler’s performance for various 
design parameters. Notice that, in the absence of channel impairments 
other than noise, the s/n is not a function of the relative powers of 
the echo and the far signal before cancellation. This result, of course, 
is not valid in a finite precision environment. 


IV. EFFECT OF PHASE JITTER AND DOUBLE-TALKING 


The distant echo, which has propagated through carrier systems, is 
likely to exhibit phase jitter. A Quadrature Amplitude Modulation 
(QAM) signal with phase jitter is usually defined as the real part of 


Z(t) = Y AnGilt — nTeloe*?), (29) 


where ©(t) is the phase jitter. 

This complex signal is generally not an analytic signal; however, for 
small enough ®&(t), it is a very good approximation to an analytic 
signal. If we assume | ®(t) | <1, we can write 


Zot) = Zi (tle = Z,(t)-[1 + jH(O)], (30) 
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where Z(t) is the analytic signal of the jitter-free, but possibly other- 
wise distorted, QAM signal. Taking the real part of (33), we get 


So(t) = si(t) — 81(t)- B(t). (31) 


It can be easily shown that s,(t) and $;(t) are uncorrelated. Therefore 
the signal s,(t)-@(t) can be considered as uncorrelated noise that is 
added to the jitter-free QAM signal s,(t). The phase jitter &(t) is 
usually slowly time varying compared to the symbol rate, so that we 
can assume that it remains constant in a given symbol period, i.e., 


(nT + iT’) = ®, = constant for given n and i = 0, i, (L — 1). 


Sampling the output signal (34) at time nT + iT’, squaring, and taking 
the average, we get 


(s§(nT + iT’)) = (si(nT + iT’)) + (s3(nT + iT’))-(®2). (32) 


If we define the average powers in a symbol period 


= 


+d (s6(nT + iT’)) (33) 


L 


1 


ne i(nT + iT’)), (34) 


oy 


nee (si(nT + iT’)) = 


Sli 
im 


Ee 
we can rewrite (36) as 
Po = P,\(1 + (#3)). (35) 


Thus the power, Po, of the far echo is generally time varying due to 
the term (®3). However, it can be shown that the power, Pi, of the 
jitter-free echo is a constant provided that L = 2. The expression, (26), 
for the residual MSE after convergence of the canceler is repeated 
here as 

ee ree C3, 

a a a TY, 
where (£’) is the average power in the uncorrelated interference signal. 
One of the components of (£’) is the quantity P,-(®2), which is due 
to the phase jitter, and another component is the desired information- 
bearing signal received from the other modem. The average power of 
this later signal, defined as in (33), will be denoted by P,. There is, of 
course, the ever-present additive Gaussian noise, and there are some 
other uncorrelated components that will generally depend on the 
specific architecture of the system. The average power of all these 
signals will be denoted by P,,. Equation (36) can now be written as 


P, + P,-(€2) + P, 
1—aNA . 


(36) 


E= (37) 
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This is the input power seen by the receiver after echo cancellation 
has taken place. The useful power for the receiver is P, and the rest 
of the power is interfering “noise”, J, defined by 


aNAP, + P;.(®2) + P, 


I=E-P,= Tn NA (88) 
The steady-state received s/n in a given symbol period is then 
P, 1—aNA 
s/n = "T= GNA + (Bi/P,)-(82) + P,P, 


Notice that the s/n is now a function of the echo-to-signal ratio 
before cancellation. The phase jitter, ®(t), is usually modeled as a 
slowly varying sine wave. As a consequence, the expression for the 
s/n is also time varying, and in order to compute the worst-case s/n, 
we have to use the largest value of (2) in (39). Assume that ®(t) is a 
simple sine wave 


@(t) = C cos 2xfot. (40) 


If B is the peak-to-peak phase jitter in degrees, then the maximum 
value of (®2) is 


2 
max (82) = (? = (3) . (41) 


This expression does not depend on frequency and can be used even 
if the jitter cannot be modeled as a sine wave. However, if we want to 
check experimental results against analytical results, we have to use 
the time average of (2), since laboratory equipment will generally 
average power measurements over long periods of time. In this case 


the quantity 
C? if¢a : 
2) — = —{ — 
(Pn) 2 2 (5 2) 2) 


should be used in (39). 


V. EFFECT OF FREQUENCY OFFSET 


The analytic signal of a QAM signal affected by frequency offset 
can be represented by 


Zo(t) = Y AnGi(t — nTheerrn, (43) 


where w, is the radian frequency offset. Frequency offset can be 
considered as a special case of phase-jitter if we replace ®(t) in (29) 
by wt. However, the expansion (30) does not hold anymore, even for 
a very small w;. In order to study the effect of frequency offset on the 
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in-band canceler we can use the analysis given in Section III and 
rewrite (43) in the following way: 


Z(t) = } ArRi(t — nT), (44) 


where, if we assume w,1' = k2z, 
A =A out (45) 


and R,(t) consists of in-phase and quadrature bandpass filters. The 
sampled output of the ith subchannel is, from (12), 


s(nT + iT’) = af’ ry; + bi V9;, (46) 


where we assume, for the time being, that &,; = 0. 
If we define 


A, = a, + jb, and Ry; = ri; + jrai, (47) 
we can rewrite (46) in the two following ways: 
2s.(nT + iT’) = (e"4A7)-R¥ + Rij. (e/"4A,)* 
= An-(e%"4Ry)* + (eo RT)-AS, (48) 


where we have defined w,T = A, and * denotes the complex conjugate. 

The two expressions in (48) lead to two different interpretations of 
the effect of frequency offset on a QAM signal. First we can consider 
that the symbol vectors A, are rotated by small increments A at the 
input of a time-invariant channel. Alternatively we can assume that 
the vectors A, are the channel’s inputs and that the sampled channel’s 
impulse responses Rj; are rotated by increments—A. The canceler’s 
performance in the presence of frequency offset is studied in Appendix 
B, where we have assumed that the nonrotated symbols A, are the 
inputs of both the canceler and the channel, and that the algorithms 
given in Section III are used for updating the tap coefficients. Under 
these conditions it is shown that the mean steady-state tap coefficients 
satisfy 

aAe!4 

1 — (1 -— aA)e*’ 
The quantity on the right is equal to the ith subchannel’s complex 
impulse response at time (n7' + 1T’) multiplied by a constant complex 
number. Therefore, the canceler’s complex taps rotate and track the 
channel’s complex impulse response with a fixed phase lag and a 
different amplitude. It is interesting to note that the mean-tap values 
go to zero when the step size a goes to zero. This most peculiar 
behavior was also observed experimentally. The ith subcanceler’s 
steady-state MSE after cancellation is derived in Appendix B. It is 
given by 


(Cri) = Rie?" (49) 
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i ay Es 1 aA_ He TR* 
(€n,i) (eni(A — 0)) oe ee aNA 072A? a (1 a qAya? euRi, (50) 
where (e2,;(A = 0)) is the MSE in the absence of frequency offset. 

The MSE for the whole canceler is obtained by averaging (50) over 
all the subcancelers. The resulting quantity is not very useful, because 
it is dependent on the channel characteristics. This is not the case for 
the Echo-Return-Loss Enhancement (ERLE), which is also studied in 
Appendix B. The ERLE is defined as the ratio between the power of 
the uncanceled echo and the power of the residual echo. It is given by 


1-—aA A? 
LE = — 9 ae eeeeecrmntcomeeans 51 
ee eee (; —aNA aA? + (1 — an) on 
This expression is the same for all the subcancelers, and therefore it 
gives also the ERLE for the whole canceler. In the derivation of (51) 
it was assumed that there was no interfering, uncorrelated signal at 
the output of the channel. 


VI. EXPERIMENTAL RESULTS 


Experimental results were obtained by using an in-band data-driven 
echo canceler implemented on a bit-slice processor. Figure 4 shows 


s/n IN DECIBELS 


@ COMPUTED 
X MEASURED 
& NONLINEARITIES 


N =32 
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STEP SIZE 


Fig. 4—The s/n in the presence of noise or nonlinearities. 
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performance curves giving the achievable s/n as a function of the step 
size for various types of channel impairments. A simulated far signal 
was added to the echo before cancellation, and the s/n was measured 
after cancellation (point A in Fig. 1b). The echo and the far signal had 
the same power before cancellation. The straight curve corresponds to 
the case where no impairments, other than linear distortion, were 
present in the echo channel. Notice that the achievable s/n, under 
these ideal conditions, increases by 3 dB when the step size decreases 
by a factor of two. One of the two curves that flatten out was obtained 
by adding noise to the echo and the far signal. In this case the noise 
power was 24 dB below the echo power. As one might expect, this 
curve goes asymptotically to an s/n of 24 dB when the step size 
becomes small. In both the preceding cases the theoretical curves were 
computed by using (28). The third curve was obtained by introducing 
nonlinearities in the echo channel. The amount of nonlinearities used 
in this experiment corresponded to the worst case reported in the 
1969-70 DDD connection survey.!° 

Figure 5 shows similar results when phase jitter was introduced in 
the echo channel. The theoretical curves were obtained by using (39) 
and (42). The frequency of the phase jitter used in the experiments 
was 120 Hz, and the far signal and far echo had again the same power. 
If the far echo’s power were X dB below the power of the far signal, 
then the flat portion of the curves would move up by X dB, as seen 


25 


20 
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Fig. 5—The s/n degradation due to phase jitter. 
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from (39). A word of caution is in order here about the usage of all the 
results discussed so far for predicting receiver performance. The noise 
in the s/n is not Gaussian, and, thus, standard probability of error 
formulas assuming Gaussian noise should not be used in a naive way. 
However, we found that the preceding results were very good indicators 
for predicting relative receiver performance for various design param- 
eters and channel impairments. One has also to consider that the 
residual echo is amplitude modulated by the phase jitter when this 
impairment is present in the echo channel, as seen from (31) and Fig. 
6. The probability of error, in this case will be mostly influenced by 
the maxima of the residual echo, rather than by its average power. 
Figures 7 and 8 show experimental results obtained when frequency 
offset was present in the echo channel. The Echo-Return-Loss En- 
hancement (ERLE) shown in Fig. 7 is defined as the ratio between 
the echo power before cancellation and the power of the residual echo 
after cancellation. No double-talker was present in these experiments. 





(b) 


Fig. 6—Residual mean-squared error (a) without phase jitter and (b) in the presence 
of phase jitter. 
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Fig. 7—ERLE degradation due to frequency offset. 
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Fig. 8—Evolution of the ERLE in the presence of frequency offset and with frozen 
taps. 


The theoretical curves were computed by using (51). The agreement 
between the experimental and analytical results is seen to be excellent. 
Notice the intuitively satisfying result that the tracking capability of 
the echo canceler improves with large step sizes. Figure 7 was obtained 
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by conducting the following experiment. The echo canceler was first 
converged when no frequency offset was present in the echo channel. 
The canceler’s taps were then frozen and a frequency offset of 2 Hz 
was introduced in the channel. Notice the periodic evolution of the 
ERLE. As was mentioned in the preceding section, the effect of 
frequency offset can be modeled by a periodic rotation of the sampled 
values of the in-phase and quadrature impulse responses of the echo 
channel. Thus, the ERLE will pass through a minimum when these 
sampled values are in phase with the complex tap coefficients of the 
echo canceler. It will pass through a maximum when the sampled 
impulse responses and the complex tap coefficients are out of phase 
by 180 degrees. 

We conclude this section with some brief comments on the impli- 
cations of the preceding results on the practical implementation of 
high-speed, full-duplex DDD modems. The effects of noise, nonlinear- 
ities, and phase jitter are very similar, as seen in Figs. 4, 5, and 7. In 
each case a comfortable s/n can be achieved by using a small enough 
step size in the tap adjustment algorithm. This is certainly true for 
the 4800 b/s modem described in Ref. 1 that requires an s/n of about 
14 dB in order to achieve a bit-error rate of 10~°. It is usually assumed 
that the far echo is always at least 10 dB below the far signal. Thus, 
under these worst-case conditions, the flat portions of the curves in 
Fig. 5 would all move up by 10 dB. In this case, even a 25-degree peak- 
to-peak phase jitter in the far echo would not deteriorate significantly 
the canceler’s performance, provided that the step size is chosen small 
enough. On the other hand, a small step size will limit the echo 
canceler’s capability of tracking frequency offset in the far echo, as 
shown in Fig. 7. A frequency offset as small as 0.01 Hz can seriously 
degrade the echo canceler’s performance. It is generally accepted that 
the United States’ domestic network does not introduce frequency 
offset in the far echo. However, this is not necessarily true in other 
countries. Frequency offset compensation techniques are known, but 
they are quite expensive to implement.*” 
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APPENDIX A 
Convergence of the MSE 


For notational convenience we will delete the index i in the equations 


but it will be understood that the analysis applies to a subcanceler and 
not the whole canceler. The updating algorithms are, from (19) and 


(20), 


Cn41 = Cn + AAnen (52) 
dizi = d, + abnen. (53) 
If we define the complex tap weights and data symbols 
C,, 2c, + jdn (54) 
A, & a, + jbn, (55) 


we can rewrite (52) and (53) in the compact form 


Chir = C, + aAnen. (56) 


Subtracting both sides of (56) from the complex vector 


R= YT + jYe, (57) 


we get the impulse-response error signal 


€nt+1 = En aAnen + F,,, (58) 


where we have defined the canceler tap errors as 


€n AR- C,, = [ri —C, + (te = d,,)]. (59) 
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In (58) we have also included a complex vector F,, which will be 
needed in Appendix B, to account for frequency offset. In this Appen- 
dix this quantity is taken to be zero. 

The canceler error at the nth iteration is, from (13), 


= (r7 rare Ci )a, + (rZ x dz )b, + L. (60) 
al TAX + ATE) + &,, (61) 


where the asterisk * denotes the complex conjugate. The MSE at the 
nth iteration becomes 


(2) =F ((eZAE + ATet)®) + (£2). (62) 


The tap error vectors e, depend only on the vectors A,-; as shown 
in (58). If we assume that consecutive data vectors An-; and A, are 
uncorrelated, and that the a;s are uncorrelated with the b/s, we can 
rewrite (62) as 


(en) = Aenex) + (En), (63) 
where A = (a2) = (b2) is the variance in the symbols. The MSE at 
the (n + 1)th iteration is 

(e741) = A(ensi€s+1) +3 (E). (64) 
Using (58), we get 
(nai) = A((en + Fi — ae, An) 


-(e& + FX — ae,A%*)) + (7) (65) 
= Alene) — aAle,(Aneg + €,A%)) 
+ 20°A?N(e2) + Af(F,), (66) 
where we have used 
ATA* = 2AN, , (67) 
and have defined 
f(Fr) = (FA(e% — aenA*) + (en — aenAn)F*). (68) 


Setting F, to zero and using (61) and (63) in (66), we get the 
equation 


(e241) = (1 — 2aA + 2a°NA?”)-(e2) + 2aA(E*). (69) 
The solution to this recurrence equation is 
(e2) = (1 — 2aA + 2a°NA”)"(e5) 
1 — (1 — 2aA + 2a°NA?’)" 


2 
i = os ORNS 0 
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APPENDIX B 
Effect of Frequency Offset 

It was shown in Section VI that the effect of frequency offset can 
be accounted for by rotating the channel’s sampled impulse responses 
at the symbol rate by angles A = w,7, where w, is the radian frequency 
offset. We can then rewrite (56) in the following way: 


Ravi — Cnai = Rn — Cyn — aAnen + Rati — Rn (71) 
€nt1 = € — AAC, + Fy, (72) 
where we have defined 
é, = R, — C, = Re’ — C, (73) 
F,, = Ras — R, = Re“™4(e“4 — 1). (74) 


The analysis given in Appendix A can now be carried through up to 
(69) by keeping F,, different from zero. The MSE at the (n+ 1)th 
iteration is now 


(@n+i) = (1 — 2aA + 2a°NA*)(en) + 2aA(é’) + Af(F,), (75) 
where from (68) 
f(Pn) = (Fie — ae, At) + (en — wen An)F%). (76) 


Using (72) and (73), we can rewrite this expression in the following 
way: 


f(Fn) = FA(RS — (Cia1)) + (RA — (Cra) ) Fa. (77) 


The evaluation of this quantity requires the knowledge of the mean- 
tap fluctuations (C,,). Replacing e, in (56) by its value in (61) and 
taking the average gives 


(Cust) = (Cn) + 5 (An(eZAS + ATe%)) 


= (C,) + aA(en). (78) 
From (73) we get 
(Chai) = (Cr)(1 — aA) + aARe~™. (79) 


After some algebra, the solution of this recurrence equation is found 
to be 


(Casi) = (Co) — aA)" 


ee —_ (1 ms aA)rteintDa 
+aA dp a 
ie 1 — (1 — aA)e”* 


In steady-state operation, as n goes to infinity, we have 


(80) 
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1 


— ep ey ee se 
(C41) aARe 1- (1 = wA)es” 


(81) 
We can now get an expression for the steady-state value of f(F,,) in 
(77). After some algebra 
f(F,) = —2R7(1 — aA) 
3 — 4 cos A + cos 2A — 2aA(1 — cos A) 


7 psa (1 — aA)? — 2(11— aA) cosA  ”’ (82) 
where we have defined 
R? = R'R*. (83) 
In the usual case where A < 1, we can use the approximation 
2 
cos A 21 - ~ (84) 
in which case (82) simplifies to 
2aAA2(1 — -R? 
f(F,) = 22S A oA) (85) 


aA? + (1 — aA)A?’ 


This quantity is a constant, which does not depend on time. The 
solution to the recurrence (75) is then, in steady-state operation, 


1-—aA AA?R? 
1 — aNA a?A? + (1 — aA)A?’ 
where (e2(A = 0)) designates the steady-state value of the MSE given 
in (70) in the absence of frequency offset. 
We can use this expression to find the value of the ERLE in the 


presence of frequency offset. We will assume that there is no other 
impairment in the channel. The output of the channel is then 


(en) = (e,(A = 0)) + (86) 


1 
Sen = 9 (A; Ri + R;,A%), (87) 


and the mean-squared output of the channel is 
(s2,) = ARZ Ri = AR’. (88) 


The MSE after cancellation is given in (86), where we put the first 
term on the right equal to zero. The ERLE is then 


ERLE = 10 log(s2,,) — 10 log(e2) 


Se a 
= 5\1 —aNA ofA? + (1— aA)A2) 


Notice that the ERLE goes to zero when the step size a goes to zero. 
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Effects of Biases on Digitally Implemented Data- 
Driven Echo Cancelers 
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(Manuscript received July 10, 1984) 


In this paper the effects of biases, of a hardware origin, on the performance 
of a digitally implemented data-driven echo canceler are studied both analyt- 
ically and experimentally. It is shown that, as a consequence of any such bias, 
the canceler tap weights can randomly drift; however, in contrast to voice- 
type cancelers and fractionally spaced equalizers, the data-driven canceler will 
not drift into instability. Nevertheless, the canceler’s performance can be 
severely degraded, even for very small amounts of bias. The main result of 
this paper is a quantitative study of the canceler’s performance as a function 
of the biases and the canceler’s various design parameters, such as the number 
of tap coefficients and the step size used in the tap adjustment algorithm. 
Although the study concentrates on the biases introduced by two’s-comple- 
ment arithmetic, the results are general enough to be used with any type of 
arithmetic, provided that the biases introduced by these different types of 
arithmetic are known. Some of the analytical results have been verified 
experimentally, in real time, on a digital signal processor constructed at AT&T 
Bell Laboratories and AT&T Information Systems. Specifically, it is shown 
how the bias introduced by rounding the product of commercially available 
two’s-complement multipliers can be eliminated by a proper choice of the 
values of the canceler’s input symbols. 


I. INTRODUCTION AND SUMMARY 


The “tap-drifting” problem in fractionally spaced equalizers and 
voice-type echo cancelers is a manifestation of the presence of small 
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biases in the digital implementation of adaptive algorithms. These 
biases allow some of the equalizer and canceler tap coefficients to 
slowly grow in magnitude. As a result, the intermediate accumulated 
sums in the filtering computations also grow. Ultimately, either the 
tap coefficients or the intermediate-accumulated sum exceeds the 
boundaries of the digital representation (“overflow”), and the algo- 
rithms become unstable. The sensitivity of these algorithms to biases 
has been investigated, and several corrective actions against tap drift- 
ing have been proposed and successfully implemented.’ Tap drifting 
can, in principle, be eliminated in its hardware origin by removing all 
the biases from the digital implementation. However, the elimination 
of biases is not always possible, especially if off-the-shelf devices are 
used in the design. For example, most commercially available arith- 
metic units use two’s-complement arithmetic, which introduces biases 
in the computations when the digital words are reduced in length. 

In this paper, we present a study of the tap-drifting problem in a 
data-driven canceler. This canceler is very attractive for two-wire full- 
duplex data transmission applications.* In its simplest configuration, 
for 4800-b/s full-duplex operation, the data-driven canceler requires 
only additions and subtractions for the computation of the filtering 
and updating algorithms. In this case no multiplier is required in the 
implementation of the canceler. However, a multiplier is required for 
higher speeds of transmission (29600 b/s). 

For most implementations, the effect of bias can be mitigated by 
using enough precision in the digital computations. However, such an 
approach will generally not be cost-effective, and in most practical 
digital implementations there will be some small biases due to the 
finite precision used in the computations. It is shown in this paper 
that, unlike voice-type cancelers,* data-driven cancelers cannot be 
driven to instability by these biases. Nevertheless, tap drifting always 
occurs to some extent and, after a sufficient period of time, introduces 
a degradation in the canceler’s performance. This phenomenon has 
been studied both analytically and experimentally. We give a general 
formula that permits the degradation in the canceler’s performance 
due to digital biases to be predicted. Real-time experimental results 
were obtained on a digital signal processor. The frequency components 
in the distortion introduced by the biases are shown to include spectral 
lines located at the origin and at multiples of the symbol rate. It is 
also shown that the performance of the biased canceler is degraded 
when the step size used in the updating of the tap coefficients de- 


* The motivations for using echo cancellation in these applications are explained in 
Ref. 3. Section II of Ref. 3 also presents more details about the echo canceler structure 
used in the following study. 
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creases. This is in direct contrast to the behavior of the unbiased, 
infinite-precision canceler for which performance improves with de- 
creasing step sizes. This phenomenon has also been observed in other 
applications.*” 

Special attention has been given to the degradation introduced by 
commercially available two’s-complement multipliers. The naive use 
of such multipliers is shown to degrade the canceler’s performance to 
unacceptable levels, even when large precision is used in the digital 
implementation. However, it is proved in this paper that proper choice 
of the values of the data symbols can completely eliminate the bias 
associated with the rounding of a two’s-complement product. Two sets 
of binary symbols having this desired property are described. 

The paper is organized as follows. In the next section the data- 
driven canceler used in the analysis and experiments is briefly de- 
scribed. In Section III we discuss the mathematical modeling of the 
biases, when two’s-complement arithmetic is utilized in the canceler’s 
updating algorithm. Quantitative results for the degradation in the 
canceler’s performance are obtained in Section IV. The frequency- 
domain characteristics of the distortion introduced by the biases are 
studied in Section V. Finally, in Section VI we present some experi- 
mental results obtained on a digital signal processor, and we compare 
these results to the analytical results. 


Il. CANCELER DESCRIPTION 


The data-driven passband echo canceler described in Refs. 8 and 3 
synthesizes a signal of the form 


s(t) = Re 3 Ang(t - net ‘ (1) 


where A,, = a, + jb, is the complex symbol to be transmitted, g(t) is a 
(possibly complex) baseband signal, w,/27 is the carrier frequency, 
and Re denotes the real part of the quantity in brackets. It is shown 
in Ref. 8 that this signal can be generated by using the structure given 
in Fig. 1. The canceler consists of two transversal filters whose taps 
are spaced at intervals T’, where 1/T” has to be at least twice the 
highest frequency in the signal s(t) in (1). This condition makes the 
canceler Nyquist, that is, it can generate an exact replica of s(t) at all 
frequencies. 

After convergence, the tap coefficients are equal to the sampled 
values of the impulse responses of the in-phase and quadrature pass- 
band filters. For this reason the canceler is called an “in-band” 
canceler, distinguishing it from other possible structures that synthe- 
size baseband equivalents of the channel. The symbol rotation at the 
input of the canceler ensures phase continuity of the carrier. For most 
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cases of practical interest, the relationship between w, and T is such 
that the rotated symbols a/, and bj, are similar to the symbols a, and 
b,. That is, if a, and b, are binary levels {+a}, then the rotated symbol 
levels will also be binary. Thus, the primes have been dropped in the 
ensuing analysis. The Digital-to-Analog (D/A) converter and low-pass 
filter at the output of the canceler perform the usual interpolation 
functions needed for further analog processing. 

Since inputs are accepted by the canceler at a rate of 1/T’, while 
the data symbols are only presented at a rate 1/T, L — 1 zero symbols, 
where L = T/T’, are inserted between successive nonzero inputs to 
the canceler. Thus, only one of every L complex taps is active for each 
filter iteration, as shown in Fig. 1. The unnecessary computations 
associated with the zero symbols can be eliminated by grouping taps 
that act simultaneously into L parallel subcancelers, as seen in Fig. 2. 
Similarly, the echo channel can be considered as a parallel combination 
of L subchannels. Convergence of the canceler is achieved by mini- 
mizing the Mean-Squared Error (MSE) between the outputs of each 
subchannel and the corresponding subcanceler. The subcancelers are 
assumed to adapt independently, and this assumption was found to be 
in excellent agreement with experimental results. The MSE for the 
whole canceler is obtained by averaging MSEs of all the subcancelers. 
A more detailed analysis of the subcanceler structure is given in Ref. 
8. The echo canceler’s performance in the presence of channel impair- 
ments is studied in a companion paper.® Some of the definitions used 
in Ref. 3 that are needed in the sequel are briefly repeated. The error 
for the ith subcanceler at time nT + iT’ is 


eni = (ri < Cn,i) Qn a (Toi ~ dni) bn + gi, (2) 


where i = 1, 2, ---, L and the superscript T denotes a transposed 
vector. In (2), ri; and ra; are the sampled in-phase and quadrature 
impulse response vectors of the ith subchannel; c,,; and d,,,; are the in- 
phase and quadrature tap coefficient vectors of the ith subcanceler; 
and ; is a signal that is uncorrelated with the data symbols (the far- 
end signal and noise). The following definitions are also needed: 


efacnt an 9 ep tt o_O 






FILTER 


bn 0 0 bn-1 0 


Fig. 1—In-band echo canceler structure. 
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SUBCANCELER (c_o,d_g) 










SUBCANCELER (c_,-1,d_)-1) 


(b) 


Fig. 2—(a) Subcanceler structure. (b) Modified canceler structure. 


An = a, + jb, = vector of complex input symbols, 


vector of complex tap coefficients 
of the ith subcanceler. 


Ga = Cn,i + Jdn,i 


The above summary of the subcanceler structure is sufficient for the 
purpose of describing the effects of biases in this paper. 


Ill. MATHEMATICAL MODEL OF BIAS 
The stochastic-gradient algorithm for the adjustment of the ith 
subcanceler’s complex tap coefficients is given by 
Ch+1i > Chai oe Ven iAny l= L 2, eiaget L, (3) 


where y is the step size, and the error e,,; is given in (2). All of the 
entries (real and imaginary) in A, are assumed to be discrete-valued 
symbols with zero mean. Biases are usually introduced in the evalua- 
tion of the correction term, 


Gai = VeniAn- (4) 
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In the following, it is assumed that this expression is computed by 
using fractional two’s-complement arithmetic, since this is the type of 
arithmetic most commonly used in digital signal processing. (Some of 
the properties of two’s-complement arithmetic are discussed in Ap- 
pendix B.) In a cost-conscious implementation, the quantity ye,,; 
should be evaluated first, because it is the same for all the tap 
coefficients. First consideration is given to the case where (4) is 
computed without utilizing a multiplier, i.e., the scaled error, ye,,;, is 
obtained by using arithmetic right shifts since y is always less than 
unity to ensure stability. Due to the finite precision of the digital 
representation, some of the lower bits of e,,; may be lost during the 
shifting operation, and consequently, a positive number will always 
decrease in magnitude and a negative number will always increase in 
magnitude. Therefore, negative bias is introduced, on the average, 
during the computation of ye,,; and during “multiplication” by data 
symbols having value less than unity (since this operation also corre- 
sponds to right shifts). If all of the symbols are binary and chosen 
equal to +1, the updating algorithm in (3) is simply implemented by 
either adding or subtracting the scaled error to each tap coefficient. 
An explicit bias occurs during these operations if the quantity ye,,; 
has “fallen out” of its register length. In this case a positive number 
becomes a true zero, but a negative number remains equal in magnitude 
to the Least-Significant Bit (LSB) of the digital representation, thus 
again introducing a negative mean bias. The correction factor in (4) 
can now be rewritten as 


Gnii = (Yeni + Aini)An + Doni = VeniAn + Bnis (5) 


where A,,,; is a random variable representing biases introduced in the 
evaluation of yep,;; Aon, is a complex random vector representing biases 
introduced in the computation of (ye,,i) An; and B,,; is the total bias. 
The other quantities are assumed to be represented with infinite 
precision. 

An expression similar to (5) is obtained when a two’s-complement 
multiplier is used to compute the correction factor in the updating of 
the tap coefficients. When a negative product is truncated, the result- 
ing number is always increased in magnitude. Conversely, when the 
product is positive, it is always decreased in magnitude. Therefore, 
truncation introduces a negative mean bias. Two’s-complement round- 
ing, on the other hand, always selects the number that is closest in 
magnitude to the exact product, independent of the product’s sign. 
One exception, however, occurs when the double-precision product is 
equally close to two single-precision numbers. In this case the most 
positive number is always selected, independently of the product’s 
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sign, and a positive bias is introduced in the computations.* This 
situation can occur quite frequently in data-driven echo cancelers 
using symbol values of +1/2 or +1/4. A remarkable property of data- 
driven cancelers is that the mean bias associated with this rounding 
can be eliminated by using the proper levels for the symbols. These 
levels must be chosen in a way such that the bias situation can never 
occur. Two sets of symbols which have the desired property are 
described in Appendix B. 

The model represented by (5) also holds for other types of arith- 
metic. Furthermore, other effects such as biases in the Analog-to- 
Digital (A/D) conversion can also be accounted for by properly defin- 
ing the random variables in Aj,,; and Ao,,;. In general, the statistics of 
these random variables will depend on 7, the statistics of e,,; and An, 
and the type of digital implementation utilized. The characterization 
of the statistics of Aj,,; and A2,,; is a formidable problem that will not 
be addressed here. However, reasonable approximations permit the 
study of the canceler’s performance degradation in the presence of 
digital biases. It will be shown later that both the mean and the 
variance of the random variables in (5) influence the canceler’s per- 
formance. However, the effect of any nonzero-mean bias will, in 
general, be predominant. This quantity is discussed next. 

From (5) the mean gradient estimate is given by 


(Gi) _ (¥en,iAn) so Bi, (6) 
where 
Bi = (Bni) = (AiniAn + Aon). (7) 


The mean values of Aj,,; and Ao,,; are generally not zero, as explained 
in the preceding discussion. In obtaining quantitative results, some 
assumptions concerning the mean values are made. First, all of the 
components of the vector (A2,,;) corresponding to the ith subcanceler 
are assumed equal. This is a reasonable assumption since the same 
scaled error is used for the updating of all of the taps of a given 
subcanceler. The weighting of this term by A, should yield the same 
average bias in steady-state operation. It is not assumed that the 
vector (Aon,;) is the same for all the subcancelers, since the error, é,,i, 
is a sample of a cyclostationary process whose statistics will depend 
upon the index 1. As a consequence, the statistics of the product, 
Yen iAn, generally vary for different subcancelers. 


* This type of bias is present in most of the commercially available two’s-complement 
multipliers. However, the bias can be removed in a new design at the cost of a slight 
increase in the chip’s complexity. This is achieved by first detecting the bias condition 
and then changing the rounding rules according to the sign of the product. 
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The influence of the term (Ai,,;A,) in (7) is somewhat more difficult 
to assess. If the assumption is made that Aj,,; is uncorrelated with A, 
and that (A,) is zero, then this quantity is also zero. However, there 
can be reasonably long sequences of symbols during which the mean 
value of A, is not zero. Biases can accumulate during these sequences 
and introduce a degradation in the canceler’s performance. We will 
not pursue this problem any further. 


IV. EXCESS ERROR DUE TO BIAS 


We will now show how biases can produce an increase in the MSE. 
It is assumed that no other impairments are present except, perhaps, 
for some uncorrelated noise added to the input signal. The mean tap 
coefficient fluctuations are investigated first. In the beginning of this 
section, the subscript i is dropped in the equations, but it will be 
understood that the analysis applied to a subcanceler. Combining (3) 
and (6), the mean tap vector evolves according to 


(Cra) = (C,) + (yenAn) + B. (8) 
Using (2), we can write this expression as 
(Cn4i) = (1 — yA)(Cn) + (m + jre)yA + B, (9) 


where A = (a2) = (b2), and the ajs and bjs are assumed to be 
uncorrelated. The steady-state tap values are given by 


B 
(Co) = tte A (10) 


where the term B/yA represents the tap deviation due to the mean 
bias B. With AC denoting the tap deviation from the optimum setting, 
one obtains 
B 
AC = (Cx) — Con = (11) 
Note that the bias is the same for each tap weight and that decreas- 
ing the step size y will result in an increased mean tap deviation. This 
contrasts with the well-known results that decreasing the step size, 
for an infinite-precision canceler (without bias), will, in general, im- 
prove steady-state performance. The deviation is also proportional to 
the bias, which agrees with intuition. Recall that in a finite-precision 
environment the step size cannot be arbitrarily close to zero, and 
therefore the tap deviation in (11) cannot approach infinity. The only 
way to effectively have a zero step size is to stop updating, in which 
case the bias no longer affects the algorithms. When quantitative 
results are discussed in the following sections, the step sizes used in 
practice will usually be found to be several orders of magnitude larger 
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than the bias B. Thus the tap offset in (11) is small. As a consequence, 
the tap coefficients will never overflow, and a data-driven canceler 
cannot become unstable. However, as will be shown later, a very small 
offset can severely degrade the MSE. 

It is interesting to contrast the preceding results with those obtained 
for a voice-type canceler, which behaves similarly to fractionally 
spaced equalizers.” In this case it can be shown that the mean tap 
offset becomes 


B 
AC = 07 =, (12) 


where . is the input autocorrelation matrix of the data signal. 
Performing a spectral expansion of the mean tap deviation yields 


AC == 5 =/fBp, (13) 
ee ie 
where y; and f; are the ith eigenvalue and eigenvector of .<% respec- 
tively. Small eigenvalues, corresponding to input frequency ranges of 
little or no energy, can make this term large, especially if {7B is not 
small. Hence, for the voice-driven canceler, distortion due to biases 
can be expected to be concentrated in frequency ranges corresponding 
to little input energy. Furthermore, the tap offset in (13) can be much 
larger than the offset given in (11), so that some tap coefficients can 
overflow and make the canceler unstable. 
The mean-squared error for the data-driven subcanceler is derived 
in Appendix A.* The expression for the ith subcanceler’s steady-state 
MSKE is given by 





N 2N 
(EB) + (82) + (Bi)? 
Y ¥ 
1—yNA , 


where (£?) is the minimum mean-squared error, N is the number of 
taps, and @; is a random variable corresponding to one component 
(real or imaginary) of the vector B,,; in (5). In the derivation of (14), 
it was assumed that (67) and (8;) were constants and the same for all 
the entries of the vector B,,;. The MSE averaged over all L subcancelers 
is 


(e2;) = (14) 


* Another expression for the MSE (also derived in Appendix A) is obtained in the 
resonance case, when the step size has its optimum value for speed of convergence. This 
step size can only be used at start-up, when there is no double-talker, in which case 
convergence is so fast that the effects of biases are negligible. Therefore, we will not 
discuss this second case any further, and the MSE discussed in the sequel corresponds 
to the nonresonant case. 
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(é) ea N (6?) ia ~ (6 ) 


(Cx) = +. (15) 
where we have defined 

y=) S er, 8 =2 5 HD (16) 

(= 7 & (8); pare 5 


and 
1 L-1 
a L » (Bi)? 


Two quantities of interest can be derived from this expression. The 
first one is the Bias-Performance Ratio (BPR), which is defined as the 
ratio of the uncanceled echo’s power to the MSE after cancellation, in 
the absence of a double-talker. This quantity is similar to the Echo- 
Return-Loss Enhancement (ERLE) studied in Ref. 3 and is particu- 
larly useful in comparing analytical and experimental results. The 
other quantity of interest is the signal-to-noise ratio (s/n), which is 
the ratio of the wanted signal power to the power of any additive noise 
after cancellation and which determines the receiver error rate. (From 
the receiver’s viewpoint, the wanted signal is the double-talker, which 
acts as noise in the tap-adjustment algorithm.) From (15) the BPR is 
given by 


P.(1 — yNA) 


LN 2LN 
= AO ag (a)? 


BPR = (17) 





where P, is the power of the uncancelled echo. In deriving an expres- 
sion for the s/n, it is assumed that the only component of the 
interfering signal is the wanted signal, i.e., there is no additive noise. 
The total signal power after cancellation is the MSE, (15), the power 
in the wanted signal is P, = (£*), and the remaining power is noise, 
so that the s/n becomes 





s/n = = = —______—___.._ (18) 
(e2) — P, EN’ og DEN. 


Under normal conditions of operation, the step size y is very small, 
relative to 1/NA, and the expressions in (17) and (18) can be approx- 
imated by 


P.Ay? 


ore © SEN(BP 


(19) 
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and 
(20) 


where it has been assumed that (8) # 0. The BPR and the s/n have 
the same expressions, except that the echo’s power, P,, in (19) is 
replaced by the wanted signal’s power, P,, in (20). Both quantities 
decrease with decreasing step size, which is in agreement with the 
preceding observations made about the mean tap behavior for small 
step size. 

Equations (14) through (20) reflect degradation in the canceler’s 
performance only under steady-state conditions. These steady-state 
quantities are not influenced by initial conditions. The complete 
equations for the MSE evolution with time are given by (55) and (64) 
in Appendix A. These equations are strongly influenced, in their 
transient terms, by the canceler’s initial state. An initial condition of 
particular interest is that which the condition that exists after the 
canceler has converged with no bias, and then a bias is introduced. As 
was discussed in Section III, this situation can arise as a result of 
certain nonrandom short-term statistics of the scrambled data se- 


quence, A,. Under these circumstances, the iterative MSE, given in 
(55), reduces to 


(é*) 
Be ee NG J 
(en) 1-— yNA 
NL 2NL 
cap (8°) se A (6)? 
a = 2A2]n 
+ (1 — [1 — 2yA + 2Ny°A‘*]”) (1 — yNA) 
4N(6*) 
ae n_ = 2A 2]n 
+ ([1 — yA]” — [1 — 2yA + 2Ny°A*)") yA(QNyA — 1)’ (21) 
and the evolution of the s/n is approximated by 
s Pyy?A 
s/t = yet ul (22) 


(e2) — P, P,(yNA) + [1 — (1 — 2yA)")(2NL)(B)?’ 


where it is assumed that the uncorrelated signal (£*) consists only of 
the wanted signal with power, P,. The approximation in (22) is 
obtained for a small step size. We will not pursue further the study of 
transient effects. 


V. FREQUENCY ANALYSIS OF BIAS DISTORTION 


For the data-driven canceler, distortion due to biases introduces 
spectral components concentrated in narrow frequency bands centered 
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around integer multiples of the symbol rate. To facilitate investigation 
- of this phenomenon, the mean tap vector is again written in terms of 
its deviation from optimum 


(Co) = Conti + AC, (23) 


where the subscript i is reintroduced to designate a subcanceler. The 
corresponding frequency response is the Fourier transform of the 
sequence obtained from the vector C.., and it is defined as 


(Cx,i(w)) = Copti(w) + ACi(w). (24) 


Since all components of the mean-bias vector B; are assumed equal 
and 


B; 
then 
N= 
AC(w) = =F eet, 
vA n=0 


3 Bi o-ieww-vy02 sinfwNT/2]} (26) 
yA sin[w7T/2] ’ 
where B; is a component of B; and 1/T is the symbol rate. 

The amplitude of the distortion is the familiar periodic sinc function. 
From (26) it is seen that the distortion is concentrated at integer 
multiples of the symbol rate. As mentioned earlier, the bias vectors 
for each subcanceler need not be the same even if the bias components 
within any particular subcanceler are assumed constant. The corre- 
sponding tap-deviation spectrum can be expressed in terms of the tap- 
deviation spectra, AC;(w), of each of the L subcancelers as 


L-1 
AC(w) = Y ACKw)e*", (27) 


where JT” = T/L. 
If each AC;(w) has the same shaping, but with a different magnitude, 
i.€., 


AC;(w) = k;ACo(w), L = 0, meas L ca 1, (28) 


where ky = 1 and the k,, are proportionality constants, then (27) 
becomes 


L-1 
AC(w) = ACo(w) YY ket". (29) 
i=0 
For the special case k, = 1 for all m, the distortion becomes 
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Fig. 3—Frequency distortion of mean tap coefficients. 


concentrated at integer multiples of 1/T’. Since all subcanceler biases 
are equal in this case, the structure becomes equivalent to one canceler 
at rate 1/7’. As the value k,, varies, bias distortion is concentrated at 
integer multiples of the symbol rate, 1/T.* One such example is shown 
in Fig. 3, where we have arbitrarily chosen values of k,, equal to 1, 0.5, 
0.25, 0.5, 1, 0.5, 0.75, and 0.5. The study of many other examples shows 
that as a general rule of thumb, the spectral lines introduced by biases 
in a data-driven echo canceler composed of L subcancelers will be 
concentrated at integer multiples of the symbol rate. Experimental 
evidence verifying these findings will be presented in the next section. 


VI. EXPERIMENTAL RESULTS 


Several experiments were conducted on a digital signal processor, 
using a 12-bit two’s-complement multiplier, to verify some of the 
qualitative and quantitative results obtained in the preceding sections. 
The echo canceler used in the experiments was operating at a sampling 
rate of 9600 samples per second, and a symbol rate of 1200 bauds. 
Thus, it could be implemented by using eight subcancelers. The 
validity of the expression for the BPR in (17) was verified by artificially 
inserting biases in the updating algorithm of an in-band echo canceler. 
A positive quantity equal to 2~ was added periodically to the correc- 
tion factor before updating the tap coefficients. Therefore, if this 
quantity was added every mT seconds, the equivalent mean bias per 


* In certain rare cases, the biases can insert nulls at some particular multiples of the 
symbol rate. 
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symbol update was (27"')/m. This bias was chosen large enough so 
that we could distinguish its effect from other digital effects such as 
round-off noise. The influence of the bias of the two’s-complement 
multiplier was eliminated by using one of the sets of symbols described 
in Appendix B. 

The measured values for the BPR are given in Fig. 4 for different 
values of m and as a function of the step size. The BPR has also been 
computed by using (17) and the corresponding curves are shown in 
Fig. 4. In Fig. 5, curves are also shown for the case in which the bias 
of the two’s-complement multiplier was not eliminated by a proper 
choice of the data symbols. No artificial bias was added, and the 
computed curves were obtained by assuming a mean bias of 27!°. (A 
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Fig. 4—Performance degradation due to simulated biases. 
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Fig. 5—Performance degradation due to rounding of a two’s-complement multiplier. 


bias of 27? was introduced whenever the bias situation occurred. With 
symbol values +1/2, this situation was likely to occur half the time so 
that a mean bias of 2-** was introduced in the algorithm.) For reason- 
ably small step sizes, both the experimental and the theoretical curves 
decrease 6 dB for each factor-of-two decrease in step size. This is 
consistent with the expression for the BPR given in (19), and similar 
behavior can be expected for the s/n as shown in (20). Both these 
quantities go to zero when the step size goes to zero. This is in direct 
contrast to the behavior of a bias-free, infinite precision canceler for 
which performance improves with decreasing step sizes. 

Although the two’s-complement multiplier rounding bias was elim- 
inated in the processor by using the sets of symbols described in 
Appendix B, there remained some other very small, unexplained 
biases. The effect of these biases on the BPR was negligible, as shown 
in Fig. 4. Nevertheless, they could be observed by studying the spec- 
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trum of the residual echo when no double-talker was present. This 
spectrum is the flat trace in Fig. 6, and the bell-shaped curve depicts 
the spectrum of the uncanceled echo. (Due to the small magnitude of 
the residual biases, it took several minutes of the canceler’s operation 
to obtain the spectrum in Fig. 6. Immediately after convergence the 
peaks were very weak.) Notice that the peaks around the origin and 
at multiples of the symbol rate (1200 bauds in this case) produce 
exactly the kind of spectrum that was predicted by the analysis in 
Section V. These peaks would be much larger if the optimum sets of 
symbols described in Appendix B were not utilized. In Fig. 6, we also 
show the spectrum of the residual echo of a voice-type canceler that 
was implemented on the processor. Notice that for the voice-canceler 
spectrum, the residual echo’s energy accumulates in the frequency 
regions where the uncanceled echo has little energy. This energy 
increases with time and ultimately, without some sort of compensation, 
the canceler diverges after several minutes of operation. Finally, for 
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Fig. 6—Power spectrum of the residual error. (a) Data-driven echo canceler. (b) 
Voice-type echo canceler. 
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Fig. 7—Bias-performance ratio for various amounts of biases. 
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Fig. 8—Signal-to-noise degradation due to biases. 


completeness, computed curves giving the BPR and the s/n for various 
values of biases are given in Figs. 7 and 8. 
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APPENDIX A 
MSE Evolution 


The following is an analysis of the evolution of the MSE for a data- 
driven echo canceler in the presence of digital biases. It is assumed 
that the input sequences a, and b, are white and uncorrelated with 
the bias, 8,, and each other. It is also assumed that the canceler spans 
the entire length of the echo impulse response. The analysis applies 
to a subcanceler, but for simplicity of notation, the index i will be 
deleted. 

The following definitions are used in the derivations: 


Bn = Bn + jBn = complex bias vector* 


C, = Cn + jdn = complex tap vector 
An = Qn + jb, = complex data vector 

+ = step size in the adjustment algorithm 

r,; = sampled in-phase channel impulse response vector (30) 


r2 = sampled quadrature channel impulse response vector (31) 


€ = sampled uncorrelated interfering signal. (32) 
The error at the nth iteration is given by 
€n = [Pr — Cn] Qn + [re — dn]7b, + &. (33) 
The in-phase and quadrature tap error vectors are defined by 
Ein Ar — Gr (34) 
Eon & To — dp (35) 


* For ease of notation the real and imaginary parts of B, are taken to be equal. 
Although this assumption is not strictly true, it will not modify the end result of the 
analysis. 
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so that (33) can be rewritten as 
Cn = ElnQn + EFnbn + &, (36) 
and the MSE becomes 
(Cnt1) = (Eins1Gne1dn+1€1n+1) 
+ (Es ntibnsibrriEenti) + (E?). (37) 


Letting 
A & (azn) = (br) (38) 


and making the usual assumption that successive data vectors are 
uncorrelated, one obtains 


(@ns1) = A(ETnt1Einti) + A(E2ntiEenti) + (E7). (39) 
The biased stochastic-gradient algorithm can be expressed as 
Cnt1 = Cn + YCndn + Bn (40) 
dni = dn + Yenbn + Bn. (41) 
Subtracting both sides of (40) and (41) from r; and r, yields 
Cit = in — Veen = Ba (42) 
EGontt = Eon — Yendn — Bn. (43) 


Using (42) and (48) in (39) gives 
(ena1) =A (Ein — YenGn — Bn) (Exn = Ve,0n = Bn)) 
+A (Eon — Yenbn — Bn)" (Ean — Yendn — Bn)) + (&?). (44) 
After some algebra, and using (36), (38), and aza, = b;b, = NA, we 
have 
(@nti) = (@n) [1 — 2yA + 2y°NA*] + 2yA(E") 
+ 2AN(Br) — 2A(Bn Ein) — 2A(BnEan). (45) 
To proceed further, (87€:) and (876) must be evaluated. We 
will assume that 6, is uncorrelated with €,,, and €2,n. 


Using (42) and (36), and assuming again that the current tap error 
vector and the current data vector are independent, we have 


(B) (Enns) = (1 — YA)(B) (En) — NCB)’, (46) 


where it is assumed that the mean vector (8,) is a constant vector 
whose entries are all equal to (8). Solving the iteration (46) yields 


1 — (1 - yA) 


(8) (Ein) = (= YA)B) Tn — SS ~N(By2, (47) 
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In (50) it is assumed that (E10) = m, that is, the canceler is originally 
started with all of its tap coefficients equal to zero. A similar expression 
holds for the quadrature component. 

The difference equation for the MSE in (45) now becomes 


(e241) = (e2)[1 — 2yA + 2Ny2A?] + 2yA(E2) + 2AN(B2) 
— 2A la — yA)"(BT)r, — 1a = 14) py 
yA 
~ 2A la — yA)“ Bn - 12 0 = 14) ip), (48) 


The above equation, although notationally complex, is only a first- 
order difference equation of the form 


(Qhsi) = (en)Ar + ki — Redd, (49) 
where 
A= 1-2yA + 2Ny7A? (50) 
hkeo=1-yA (51) 
hy = 2yA(E2) + 2AN(B?) + ~N(BY* (52) 
ho = 2A(BT)r, + 2A(BT)rp + ~ (BY (53) 


and where it is assumed that (62) is a constant and equal to ({7). 
Assuming stability, 


JAif <1 and |dA:| <1, (54) 


there are two cases to be considered in solving (53). 
Case 1: \; ¥ As (nonresonant). 
Assume 


(e2) => C, + Covt + C33 ‘. (55) 


Then solving (49) by substitution, and assuming the initial condition 
(e%) gives 











C=, (56) 
_ ke ke 

ae Pe Cia Te en 

Cy = (e6) = C; <— Cs . (58) 


In terms of the original parameters in (48), 
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(2) +2 (6%) + 2 ay? 
‘y. 


= Y 
C, = = ae aye > NA) (59) 
4N 
= (8) n+ (Brad + NB)? 
C; = eres Oe _ ee (60) 
(2NyA — 1) 
Cr = (e8) - G.-C. (61) 


All of the parameters appearing in the expression for the MSE in (55) 
are now known. In the steady state, as n goes to infinity the MSE 
becomes 


N 2N 
(7) + a (8?) + 7A (By? 


(e2) = 1-yNA (62) 
Case 2: Resonance (71 = 72). 
In this case, from (50) and (51) 
= NT (63) 
Assume 
(eh) = C, + Coy? + nCeyt. (64) 


Solving (49) again by substitution, we get, after some algebra, 
(en) = 2(£") + 4NA(B”) + 16N°A(B)? 


+ [(e3) — 26°) — 4NA(B?) — 16N*A(B)”] ( = ax) 


— n[2A[(B")r + (8")r + 8N7A (6) (1-4) “. 5) 
In this case the steady-state MSE is 
(e2) = 2(£°) + 4NA(B?) + 16N?7A(B)?. (66) 


APPENDIX B 
Bias-Free Two’s-Complement Multiplication for a Data-Driven Canceler 


In this Appendix, it is shown that the bias that occurs when rounding 
the product of a two’s-complement multiplier can be eliminated by 
properly choosing the values of the binary input symbols. Let a be an 
(/+1)-bit fractional two’s-complement number represented as 


Q = ApMAz +++ A, (67) 
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where a; = 0 or 1, and qa; is the LSB of the number. 
The numerical value of a is given by* 


l 
a=-aA + » a;2~*, (68) 


i=1 


and the product of two such numbers becomes 


I I 
ax = (—ay + ¥ a;27')(—xo + YS x,27*) 
Zt k=l 


l l Lol 
= ApXp — Ay YD x2" — x Ya2i+ ¥% Y ax,2-“”. (69) 
k=1 i=1 i=1 R=1 
This can be rewritten as a two’s-complement number 

al 


ax = to + ¥ ¢274, (70) 
i=1 

where the summation on the right now goes to 2/, making the length 
of the product (2/ + 1). A two’s-complement multiplier usually has one 
of two means of reducing this product to a (J + 1)-bit number. In 
truncation, all the bits corresponding to i => /+ 1 are discarded. From 
(70), it is seen that such an operation always decreases the magnitude 
of a positive number, and increases the magnitude of a negative 
number, thus introducing a negative bias. In rounding a number, 27? 
is added to the product in (70) and the result is truncated to (J + 1) 
bits. This always selects the (J + 1)-bit number that is closest in 
magnitude to the true product. An ambiguity arises, however, when 
this product is equidistant from two (J + 1)-bit numbers. It is seen 
from (69) that rounding, in this case, always increases the magnitude 
of a positive number and decreases the magnitude of a negative 
number, thus introducing a positive bias in the arithmetic. 

It is assumed that the numbers +a and x in (69) represent the 
symbols and the scaled error, respectively, in the updating algorithm. 
The bias situation in rounding arises when the following conditions 
occur in (70): 


C41 =1 and ¢;=0 for i>/l+1. (71) 
Inspection of (69) shows that these conditions are equivalent to 
I l 
DS ae, 2 = gen, (72) 
i=1 k=U+1-i 


* It is readily verified that, with this definition, the largest negative number that can 
be represented is —1, and the largest positive number is +1 — LSB. 
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The rounding bias can be eliminated by choosing the a,’s in such a 
manner as to never satisfy (72), regardless of the choice of the x;,’s. 
The search for all of the possible numbers, a, having this property is 
tedious and will not be pursued here. Rather, two sets of symbols are 
proposed and shown to have the desired property. The proof requires - 
that the scaled error be reasonably small. More specifically, the mag- 
nitude of x must be strictly less than one-half. This assumption is 
always satisfied in practice, even during start-up. 
The first two binary symbols to be considered are defined by 


a = 0100 --- 001 = 27+ 27 (73) 


and 


l 
-a=1011---1l1l=-1+ ¥ 27% (74) 
i=2 
These numbers are equal to {+((1/2) + LSB)}. Considering the positive 
symbol and replacing the a,’s by their value in (72) yields 


l 
aj gO) = OD, (75) 
k=1 


This equation can only be satisfied under the condition x, = 1 and 
x; = 0 for i > 1. From (68), it is seen that the only two numbers 
satisfying this condition are +1/2. These numbers were discarded 
earlier as valid solutions. For a negative symbol (72) becomes 


l l 
y SS aoe Soe. (76) 
i=2 k=l+1-i 
A sequence of x;’s cannot be synthesized to satisfy this equation. The 
lowest-order bit on the left corresponds to the power 2I of 1/2. There 
is only one term on the left in (76) contributing to this bit’s value: 1 = 
k = I. Since there is no such term on the right, this bit must be zero 
(x; = 0). The value of second-lowest bit corresponding to a power 
(2I-1) depends upon two terms: i = 1,k = 1—1andi=!-—1. Therefore, 
its value is (x; + x-1), which must be zero modulo-2. Since x; = 0, this 
implies x,_; = 0. Reasoning by induction shows that all of the x,;’s must 
be zero for i = 2. The highest-order bit corresponds to the (J + 1)th 
power of 1/2. Combining this with the preceding result, we see that 
only one term contributes to its value; 1 = 1, k = 1. Therefore, to satisfy 
(76), x; must equal one. As previously noted, the only two numbers 
that produce the bias situation are 1/2, which are not valid solutions. 
The second set of symbols that eliminates the multiplier’s bias is 
given by 


a=0O111---111=1-—-27 (77) 
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and 
—a = 1000 --- 001 = —1 + 2! (78) 


These numbers are equal to {+(1-LSB)}. Proof by contradiction, 
similar to the one utilized previously, can be used to show that these 
symbols eliminate the multiplier’s bias. The proof is left as an exercise 
for the reader. 


AUTHORS 


John M. Cioffi, B.S. (Electrical Engineering), 1978, University of Illinois; 
M.S. and Ph.D. (Electrical Engineering), Stanford University, in 1979 and 
1984, respectively; Bell Laboratories, 1979-81; IBM Research, 1984—. In 1978 
Mr. Cioffi received the Jordan and Davidson awards from the University of 
Illinois as the outstanding senior in Engineering and the Valedictorian of 
Engineering, respectively. He is the author of numerous publications in the 
areas of adaptive filtering and data transmission and storage. Member, Tau 
Beta Pi, Eta Kappa Nu, Phi Kappa Phi, Phi Eta Sigma, Sigma Xi, IEEE. 


Jean-Jacques Werner, Ing. Deg., 1965, INSA, Lyon, France; M.S. (Electri- 
cal Engineering), 1967, Laval University, Canada; Eng. Sc.D. (Electrical 
Engineering), 1973, Columbia University; Bell Laboratories, 1973-1982; 
AT&T Information Systems, 1983—. At Bell Laboratories and AT&T Infor- 
mation Systems, Mr. Werner has worked on problems in data transmission 
and digital signal processing. Member, Sigma Xi, IEEE. 


t 


138 TECHNICAL JOURNAL, JANUARY 1985 


PAPERS BY AT&T BELL LABORATORIES AUTHORS 


COMPUTING/MATHEMATICS 


Chang J. M., Maxemchuk N. F., Reliable Broadcast Protocols. ACM T Comp 
2(3):251-273, Aug 1984. 

Cody W. J. et al., A Proposed Radix-Independent and Word-Length-Independent 
Standard for Floating-Point Arithmetic. IEEE Micro 4(4):86-100, Aug 1984. 
Fishburn P. C., Elements of Risk Analysis in Non-Linear Utility Theory. Infor 
22(2):81-97, May 1984. : 

Kapur R. N., Browne J. C., Techniques for Solving Block Tridiagonal Systems 
on Reconfigurable Array Computers. SIAM J Sci 5(3):701-719, Sep 1984. 

Pierce J. R., More on the Fifth Generation (Letter). Abacus-NY 1(4):8, Sum 1984. 
Prell KE. M., Sheng A. P., Building Quality and Productivity Into a Large Software 
System. IEEE Softw 1(3):47-54, Jul 1984. 

Whitt W., The Amount of Overtaking in a Network of Queues. Networks 
14(3):411-426, Fal 1984. 


ENGINEERING 


Acampora A. S., Hluchyj M. G., A New Local Area Network Architecture Using 
a Centralized Bus. IEEE Comm M 22(8):12-21, Aug 1984. 

Alferness R. C., Buhl L. L., Low-Crosstalk Waveguide Polarization Multiplexer/ 
Demultiplexer for \ = 1.32 um. Optics Lett 9(4):140-142, 1984. 

Arnold H. W., Bodtmann W. F., Switched-Diversity FSK in Frequency-Selective 
Rayleigh Fading. IEEE Veh T 33(3):156-163, Aug 1984. 

Baumert R. J., Cameron L. E., Wilson R. A., A Mixed EFL I?L Digital Telecom- 
munication Integrated Circuit. IEEE Device 31(2):160-165, 1984. 

Chin A. K., Caruso R., Young M.S. S., Vonneida A. R., Uniformity Characterization 
of Semi-Insulating GaAs by Cathodoluminescence Imaging. Appl Phys L 
45(5):552-554, Sep 1, 1984. 

Chraplyvy A. R., Marcuse D., Henry P. S., Carrier-Induced Phase Noise in Angle- 
Modulated Optical-Fiber Systems. J Lightw T 2(1):6-10, Feb 1984. 

Fisher R. E., UHF Television Interference Associated With Cellular Mobile 
Telephone Systems. IEEE Veh T 33(3):244-249, Aug 1984. 

Forrest S. R., Gain-Bandwidth-Limited Response in Long-Wavelength Ava- 
lanche Photodiodes. J Lightw T 2(1):34-39, Feb 1984. 

Haight R., Bokor J., Storz R. H., Stark J. B., Bucksbaum P. H., Freeman R. R., 
Photoemission Apparatus Using XUV Harmonics of a Picosecond KRF 
Laser. P Soc Photo 476:61-64, 1984. 

Hartman D. H., Multiple Reflection Noise Production on Transmission-Systems 
With Periodic or Quasi-Periodic Impedance Disturbances—-A Time-Domain 
Approach. IEEE Circ S 31(10):866-875, Oct 1984. 

Hayes J. R., Leheny R. F., Temkin H., Gossard A. C., Wiegmann W., Electrolumi- 
nescence From a Heterojunction Bipolar Transistor. Appl Phys L 45(5):537- 
539, Sep 1, 1984. 

Levy U., Logan R. A., Niv Y., Laser Cathode-Ray Tube Operation at Room 
Temperature. Appl Phys L 45(5):497-499, Sep 1, 1984 . 

Linke R. A., Direct Gigabit Modulation of Injection Lasers—Structure-De- 
pendent Speed Limitations. J Lightw T 2(1):40-43, Feb 1984. 

Liu P. L., Ogawa K., Statistical Measurements as a Way to Study Mode Partition 
in Injection Lasers. J Lightw T 2(1):44-48, Feb 1984. 

McCaughan L., Low-Loss Polarization-Independent Electrooptical Switches at 
Lambda = 1.3 um. J Lightw T 2(1):51-55, Feb 1984. 


Seth S. C., Agrawal V. D., Characterizing the LSI Yield Equation From Wafer 
Test Data. IEEE Comp A 3(2):123-126, Apr 1984. 


139 


Stone J., Marcuse D., Direct Measurement of Second-Order Dispersion in Short 
Optical Fibers Using White-Light Interferometry. Electr Lett 20(18):751-752, 
Aug 30, 1984. 

Winters J. H., Optimum Combining in Digital Mobile Radio With Cochannel 
Interference. IEEE Veh T 33(3):144-155, Aug 1984. 

Yeh Y. S., Wilson J. C., Schwartz S. C., Qutage Probability in Mobile Telephony 
With Directive Antennas and Macrodiversity. IEEE Veh T 33(3):123-127, Aug 
1984. 


MANAGEMENT/ECONOMICS 


Doshi B. T., Lipper E. H., The Throughput Performance of a Prioritized LIFO 
Service Discipline. Oper Res L 3(2):75-80, Jun 1984. 


PHYSICAL SCIENCES 


Batlogg B., Remeika J. P., Cooper A. S., Fisk Z., Magnetism and Superconductivity 
in CeCu.Siz Single Crystals. J Appl Phys 55(6):2001-2003, 1984. 

Belfiore L. A., Schilling F. C., Tonelli A. E., Lovinger A. J., Bovey F. A., Magic Angle 
Spinning Carbon-13 NMR Spectroscopy of Three Crystalline Forms of Poly- 
butene-1. Polym Prepr 25(1):351-353, 1984. 

Beni G., Hackwood S., Intermittent Turbulence and Period Doubling at the 
Corrosion Passivity Transition in Iron. J Appl Elec 14(5):623-626, Sep 1984. 


Benton J. L., Levinson M., Macrander A. T., Temkin H., Kimerling L. C., Recombi- 
nation Enhanced Defect Annealing in N-InP. Appl Phys L 45(5):566-568, Sep 1, 
1984 . 

Bergmann E. E., McCaughan L., Watson J. E., Coupling of Intersecting Ti-LiNbO3 
Diffused Waveguides. Appl Optics 23(17):3000-3003, Sep 1, 1984 . 

Caton J. A., Heywood J. B., Mendillo J. V., Hydrocarbon Oxidation in a Spark- 
Ignition Engine Exhaust Port. Comb Sci T 37(3-4):153-169, 1984. 


Davis G. P., Moore C. A., Gottscho R. A.. Dynamics of Laser Stimulated Etching 
of Germanium by Bromine. J App! Phys 56(6):1808-1811, Sep 15, 1984. 


Donnelly V. M., Karlicek R. F., Excimer Laser Enhancement and Probing of III- 
V Compound Semiconductor Chemical Vapor Deposition. P Soc Photo 476:102- 
109, 1984. 


Dutta N. K., Napholtz S. G., Yen R., Brown R. L., Shen T. M., Olsson N. A., 
Craft D. C., 1.3 um InGaAsP DCPBH Multiquantum-Well Lasers. Electr Lett 
20(18):727-728, Aug 30, 1984. 

Freund R. S., Donohue D. E., Fisanick G. J., Polarization Line Shape of Balmer-8 
From Electron Impact Dissociation of H?.. J Chem Phys 80(5):1754-1759, 1984. 
Gershenfeld N., The Measurement of Optically Thick Atomic Vapor Densities 
by the Nonlinear Least-Squares Fitting of Absorption or Fluorescence Spectra 
(Letter). Nucl Inst A 224(3):570-572, Jul 15, 1984 . 

Harrus A., Mihalisin T., Batlogg B., Saturation of Resistivity in the La;.,Ce,Rhe 
Mixed Valent Kondo System. J Appl Phys 55(6):1993-1995, 1984. 

Harrus A., Timlin J., Mihalisin T., Batlogg B., Resistivity, Susceptibility, and 
Superconductivity in Mixed Valent Ce(Rh,.,Ru,)2”. J Appl Phys 55(6):1990- 
1992, 1984. 

Hodeau J. L., Marezio M., Remeika J. P., The Structure of [Er(1);.,, Sn(1)x] 
Er(2)4RhgSn(2)48n(3)12Sn(4)2, a Ternary Reentrant Superconductor. Act 
Cryst B 40(Feb):26-38, 1984. 

Johnson A. M., Glass A. M., Olson D. H., Simpson W. M., Harbison J. P., High 
Quantum Efficiency A-Si-H Picosecond Transit Time Limited Schottky-Bar- 
rier Photodetectors. J Non-Cryst 66(1-2):381-386, Jul 1984. 

Kahane C., Frerking M. A., Langer W. D., Encrenaz P., Lucas R., Measurement of 
the Formaldehyde Ortho to Para Ratio in Three Molecular Clouds. Astron 
Astr 137(2):211-222, Aug 1984. 


140 TECHNICAL JOURNAL, JANUARY 1985 


Kash K., Shah J., Carrier Energy Relaxation in Ino.53Gao,47 as Determined From 

Picosecond Luminescence Studies. Appl Phys L 45(4):401-403, Aug 15, 1984 . 

Lang D. V., Cohen J. D., Harbison J. P., Chen M. C., Sergent A. M., Resolution of 

the A-Si-H DLTS Energy Scale Controversy. J Non-Cryst 66(1—2):217-222, Jul 

1984. 

Levy R. A., Green M. L., Gallagher P. K., Characterization of LPCVD Aluminum 

for VLSI Processing. J Elchem So 131(9):2175-2182, Sep 1984. 

Mazumdar P., Bhagat S. M., Manheimer M. A., Chen H. S., Low Field Magnetization 

Studies on (Fe,Nij-x)7sPieBeAls Ribbons. J App! Phys 55(6):1685-1687, 1984. 

McNevin S. C., Becker G. E., CF4/Silicon Surface Reactions: Evidence for Parallel 

Etching Mechanisms From Modulated Ion Beam Studies. J Vac Sci B 2(1):27- 

33, 1984. 

Mims W. B., Elimination of the Dead-Time Artifact in Electron Spin-Echo 

Envelope Spectra. J Magn Res 59(2):291-306, Sep 1984. 

Murarka S. P., Effect of Oxygen Contamination on the Properties of Cosputtered 

Tantalum Silicide. Appl Phys L 45(4):392-394, Aug 15, 1984. 

Smith J. L., Fisk Z., Willis J. O., Batlogg B., Ott H. R., Impurities in the Heavy- 

Fermion Superconductor UBe,3. J Appl Phys 55(6):1996-2000, 1984. 

Sodini C. G., Ko P. K., Moll J. L., The Effect of High Fields on MOS Device and 

Circuit Performance. IEEE Device 31(10):1886-1393, Oct 1984. 

Stern M. B., Craighead H. G., Liad P. F., Mankiewich P. M., Fabrication of 20-nm 

Structures in GaAs. Appl Phys L 45(4):410-412, Aug 15, 1984 . 

Tonelli A. E., Schilling F. C., C-NMR Chemical Shifts and the Microstructure 

of Propylene-Vinyl Chloride Copolymer With Low Propylene Content. Polym 

Prepr 25(1):334-335, 1984. 

Tyson J. A., Comparison of Space Telescope and 4-Meter Ground-Based Tele- 

ee ane Galaxy Detection and Photometry. Pub Ast S P 96(581):566-573, 
ul 1984. 

Weschler C. J., Indoor Outdoor Relationships for Nonpolar Organic Constituents 

of Aerosol Particles. Env Sci Tec 18(9):648-652, Sep 1984. 

Weschler C. J., Sulfur-Dioxide Content of Mount St. Helens Ash. J Geo Res-A 

89(ND3):4891-4894, Jun 20, 1984 . 

Wicksted J. P., Shapiro S. M., Chen H. S., Investigation of the Ferromagnetic- 

Spin Glass Transition in a-(Fe77Cre3)75PiesBeAls. J Appl Phys 55(6):1697-1699, 

1984. 

Yamada Y., Fujii Y.. Akahama Y., Endo S., Narita S., Axe J. D., McWhan D. B., 

Lattice-Dynamical Properties of Black Phosphorus Under Pressure Studied 

by Inelastic Neutron Scattering. Phys Rev B 30(5):2410-2413, Sep 1, 1984. 

Zucker J. E., Pinczuk A., Chemla D. S., Gossard A., Wiegmann W., Optical Vibrational 

Modes and Electron-Phonon Interaction in GaAs Quantum Wells. Phys Rev L 

53(13):1280-1283, Sep 24, 1984. 


SOCIAL AND LIFE SCIENCES 


Shute S. J., Starr S. J., Effects of Adjustable Furniture on VDT Users. Human 
Fact 26(2):157-170, Apr 1984. 


PAPERS BY AT&T BELL LABORATORIES AUTHORS 141 


CONTENTS, FEBRUARY 1985 


Part 1 


Algorithms for Estimation of Three-Dimensional Motion 
J. Salz and A. N. Netravali 


Homenet: A Broadband Voice/Data/Video Network on CATV 
Systems 
M. Hatamian and E. G. Bowen 


Analysis of a Multistage Queue 
B. T. Doshi and K. M. Rege 


A Probabilistic Distance Measure for Hidden Markov Models 
B.-H. Juang and L. R. Rabiner 


A Conditional Response Time of the M/M/1 Processor-Sharing Queue 
B. Sengupta and D. L. Jagerman 


A Study on the Ability to Automatically Recognize Telephone-Quality 
Speech From Large Customer Populations 
J. G. Wilpon 


Part 2 
COMPUTING SCIENCE AND SYSTEMS 


Unintrusive Communication of Status in a Packet Network in Heavy 
Traffic 
G. J. Foschini 


A Note on Discrete Representation of Lines 
M. D. Mellroy 


Models for Configuring Large-Scale Distributed Computing Systems 
B. Gavish 


Inverted Decision Tables and Their Application: Automating the 
Translation of Specifications to Programs 
L. S. Levy and H. T. Stump 


Proposed Specification of BX.25 Link Layer Protocol 
R. P. Kurshan 


Approximate Analysis of a Generalized Clocked Schedule 
A. A. Fredericks, B. L. Farrell, and D. F. DeMaio 


143 


Numerical Computation of Delays in Clocked Schedules 
M. H. Ackroyd 


Analysis of Clocked Schedules—High-Priority Tasks 
B. T. Doshi 


144 TECHNICAL JOURNAL, JANUARY 1985 





AT&T TECHNICAL JOURNAL is abstracted or indexed by Abstract Journal in 
Earthquake Engineering, Applied Mechanics Review, Applied Science & Technology 
Index, Chemical Abstracts, Computer Abstracts, Current Contents/Engineering, 
Technology & Applied Sciences, Current Index to Statistics, Current Papers in Electrical 
& Electronic Engineering, Current Papers on Computers & Control, Electronics & 
Communications Abstracts Journal, The Engineering Index, International Aerospace 
Abstracts, Journal of Current Laser Abstracts, Language and Language Behavior 
Abstracts, Mathematical Reviews, Science Abstracts (Series A, Physics Abstracts; Series 
B, Electrical and Electronic Abstracts; and Series C, Computer & Control Abstracts), 
Science Citation Index, Sociological Abstracts, Social Welfare, Social Planning and 
Social Development, and Solid State Abstracts Journal. Reproductions of the Journal 
by years are available in microform from University Microfilms, 300 N. Zeeb Road, 
Ann Arbor, Michigan 48106. 


rN 


AVS 


